Towards an Effective Spin Hamiltonian of the Pyrochlore Spin Liquid Tb2Ti207 
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Tb2Ti2 07 is a pyrochlore antiferromagnet that has dynamical spins and only short-range corre- 
lations even at 50 mK — the lowest temperature explored so far — which is much smaller than the 
scale set by the Curie- Weiss temperature ^cw ~ — 14 K. The absence of long-range order in this 
material is not understood. Recently, virtual crystal field excitations (VCFEs) have been shown to 
be significant in Tb2Ti2 07, but their effect on spin correlations has not been fully explored. Building 
on the work in Phys. Rev. Lett. 98, 157204 (2007), we present details of an effective Hamiltonian 
that takes into account VCFEs. Previous work found that VCFEs-induced renormalization of the 
nearest neighbor Ising exchange leads to spin ice correlations on a single tetrahedron. In this paper, 
we construct an effective spin-1/2 low-energy theory for Tb2Ti207 on the pyrochlore lattice. We 
determine semiclassical ground states on a lattice that allow us to see how the physics of spin ice 
is connected to the possible physics of Tb2Ti2 07. We observe a shift in the phase boundaries with 
respect to those of the dipolar spin ice model as the quantum corrections become more significant. 
In addition to the familiar classical dipolar spin ice model phases, we see a stabilization of a q = 
ordered ice phase over a large part of the phase diagram — ferromagnetic correlations being pre- 
ferred by quantum corrections in spite of an antiferromagnetic nearest neighbor exchange in the 
microscopic model. Frustration is hence seen to arise from virtual crystal field excitations over and 
above the effect of dipolar interactions in spin ice in inducing ice-like correlations. Our findings 
imply, more generally, that quantum effects could be significant in any material related to spin ices 
with a crystal field gap of order 100 K or smaller. 

PACS numbers: 75.10.Dg, 75.10.Jm, 75.40.Cx, 75.40.Gb 



I. INTRODUCTION 



The problem of finding a low energy effective theory 
from a microscopic theory or directly from experimental 
considerations is a ubiquitous one in physics. The pur- 
pose is to identify the relevant degrees of freedom at some 
energy scale in order to capture the important physics at 
that scale. Often in condensed matter physics, a large 
separation of energy scales facilitates the process of find- 
ing an effective theory: for example in the spin icesi*^!^ 
discussed below. When the separation of scales is not 
large, virtual (quantum mechanical) processes can be- 
come important, as in the Kondo problem in which dou- 
ble occupancy of the impurity in the Anderson model can 
be treated as a virtual process that generates the well- 
known s-d exchange interaction4. One focus of this paper 
is the construction of such a low energy effective theory 
for a highly exotic magnetic material - the Tb2Ti207 py- 
rochlore magnetic material. 

A second thread to the present work is frustration, 
which occurs in magnetism when interactions between 
spins cannot be minimized simultaneously. This hap- 
pens, in the case of geometric frustration, as a conse- 
quence of the topology of the lattice. As an exam- 
ple, antiferromagnetic isotropic exchange interactions be- 
tween classical spins on the vertices of the three di- 
mensional pyrochlore lattice of corner-sharing tetrahedra 
are frustrated. ^^'^'^'^'^ One consequence of this frustra- 
tion is an extensive (macroscopic) ground state degen- 
eracy and lack of conventional long-range order down 



to arbitrarily low temperatures. Theoretically, this de- 
generacy is expected to be lifted, partially, or fully, by 
other interactionsr^ii^ perhaps assisted by the presence 
of thermal or quantum fluctuations il d'^d^ These lessons 
carry over to real pyrochlore magnets in which the frus- 
tration of the principal spin-spin interaction usually man- 
ifests itself in a transition to long-range orde r^^d^d^d^ 
or a spin glass transitio n^^d^ well below the temperature 
scale set by the interactions — the Curie- Weiss temper- 
ature In fact, this is a ubiquitous fingerprint of 
highly frustrated magnets. 

When short-range spin correlations persist down to ar- 
bitrarily low temperatures, as in the isotropic exchange 
pyrochlore antiferromagnet of Refs. 00, the system is 
referred to as a spin liquid or collective paramagneti^ 
Given the large proportion of geometrically frustrated 
magnetic materials which have been studied experimen- 
tally and which do ultimately exhibit an ordering transi- 
tion, it does seem that spin liquids are rather rare in 
two and three dimensions j ^^i^^'^^i^^'^'^i^^'^^ One would 
expect, on general grounds, this scarcity to be partic- 
ularly apparent in three dimensional materials where 
thermal and quantum fluctuations are the most eas- 
ily quenched. This paper is concerned with the mate- 
rial Tb2Ti207 which is one of the very few three di- 
mensional spin liquid candidates!^ Tb2Ti207 is a py- 
rochlore antiferromagnet that is not magnetically ordered 
at any temperature above the lowest explored tempera- 
ture of 50 mK,— although the Curie- Weiss tem- 
perature, 6'cw, is about —14 K, that is three hundred 
times larger)^ Despite ten years^^ of experimental and 
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FIG. 1: (color online). Cubic unit cell of the pyrochlore lat- 
tice. The spin configuration shown is the ordered LRSIooi 
state of the dipolar spin ice model.— The spins on each tetra- 
hedron are aligned in the local [Iff] direction and satisfy the 
two-in/two-out ice rule. 



theoretical interest in this system, the low energy mag- 
netic properties of this material are still not currently 
understood4ii^i2^ 

In this article, we build on earlier work^ by present- 
ing further evidence that qualitatively new physics, in 
the form of geometrical frustration, is generated via vir- 
tual crystal field excitations (VCFEs) in TbsTisOr. The 
frustration of interactions coming from high energies is 
not without precedent in condensed matter physics: frus- 
trated exchange beyond nearest neighbor and ring ex- 
change terms arise in small t/U effective theories derived 
from the Hubbard model at half-fillingi^i^^i^ In this 
problem, the higher order terms in the effective model 
have only a quantitative effect on the physics which is 
already captured by the lowest order terms.— 

In contrast, qualitatively new phenomena have been 
proposed to arise by integrating out high energies in a 
recent work on Mott systems?^ and in gauge theories 
of frustrated magnetic systems^i^ (which, interestingly, 
take as starting points models closely related to the ef- 
fective model derived in Sections IIIII and IIVI of this pa- 
per). The substantial effect of VCFEs on low energy 
physics advocated in Ref . [s^ and in this article is reminis- 
cent of the recent experimentally motivated proposal that 
PrAu2Si2 is a disorder- free spin glass owing to frustration 
dynamically arising from excited crystal field levels.— 
Before launching into the calculations, we first describe 
some earlier developments relating to Tb2Ti207 to mo- 
tivate our approach to this problem. 



A. Phenomenology of Tb2Ti207 

There is one particular property that may be useful for 
making progress towards understanding the low energy 
physics of Tb2Ti207 and which is shared by all the com- 
pounds in the R2M2O7 family of compounds to varying 
degreesi^ (here B?^ is a rare earth ion with a magnetic 
crystal field ground state and M'^+ is non-magnetic Ti''+ 
or Sn'*+). It is the smallness of the energy scale due to 
interactions, V , compared with the crystal field splitting, 
A, between the single ion ground state doublet and the 
first (lowest) excited states. The interactions are typi- 
cally of the order of 0.1 K or smaller while the lowest 
crystal field splitting is of the order of tens or hundreds 
of Kelvini^i^i^ This means that the ground state wave- 
function and low energy excitations mainly "live" in the 
Hilbert space spanned by the ground state crystal field 
states on all lattice sites. As we shall see in detail later 
on, the interactions, V, admix excited crystal field wave- 
functions into the ground state doublet and these quan- 
tum corrections are weighted by (y)/A4^ For the spin 
ices, IIo2Ti207 and Dy2Ti207, for which A is of the or- 
der of 300 the effect of excited crystal field levels 
can be ignored to a very good approximation and the 
angular momenta can then be treated as classical Ising 
spins jii^i^ In common with the spin ices, Tb2Ti207 has 
a crystal field ground state that can be described in terms 
of Ising spins But, the (classical) dipolar spin ice model 
(DSIM) which has, through various studies demonstrated 
its veracity in comparisons to the spin ices^^i^i^ is not 
a good model for Tb2Ti207. 

An estimate of the antiferromagnetic exchange cou- 
pling in Tb2Ti207 — puts this compound close to the 
phase boundary of the DSIM between the paramagnetic 
spin ice state (or lower temperature long-range ordered 
spin ice phase) and the four sublattice long-range Neel 
antiferromagnetic phase (see inset to Fig. [5])'^'^'^ None 
of these states adequately describes Tb2Ti207. The long- 
ranged ordered phases can be ruled out on the grounds 
that no Bragg peaks are observed in the diffuse neutron 
scattering patternj ^^i^^ A comparison with spin ice phe- 
nomenology is a little more subtle. One of the main fea- 
tures of the spin ice state is that it harbors a large resid- 
ual entropy as deduced by integrating the heat capacity 
downwards from high temperatures 4^ Whereas, similarly 
to what has been observed in spin ices^''^ there is a broad 
bump in the specific heat Cy between 1 K and 2 K as 
the temperature is lowered, at present it remains diffi- 
cult to determine whether there is a residual entropy in 
the collective paramagnetic state of Tb2Ti207 j'^"i^" The 
study in Ref. |50 finds a slightly different heat capacity 
to the one in Ref. [13 and claims no evidence of residual 
entropy in Tb2Ti207 owing to almost a complete recov- 
ery of the full entropy of the doublet-doublet crystal field 
levels (see also Ref. [5l| for a similar finding). Instead it 
reports that there is a sharp feature in the heat capacity 
at about 300 mK indicating the onset of a glassy state. 
Classiness has also been observed in the susceptibility 



3 



0.06 



0.05 



0.04 



•0.03 



<l 



0.02 



0.01 



LRSI 




PM 




Spin Ice 

1 


AIAO 


LRSI 001 





-0.3 -0.2 -0.1 0.1 0.2 0.3 

Xx (,K) 



0.1 0.2 0.3 

Jex (K) 



0.4 



FIG. 2: (color online). Semiclassical ground state phases for 
the cubic unit cell model with Ewald summed dipole-dipole 
interactions as the crystal field gap, A, and the bare exchange 
coupling, e7cx, are varied. The horizontal bar indicates a value 
for 1/A (A = 18 K) and a range of ^/ex that are consistent 
with experimental results on Tb2Ti207.— The inset is the 
phase diagram of the dipolar spin ice model^i^ adopted for 
Tb2Ti207 with V = 0.0315 K with a vertical dotted line 
showing an estimated ^cx = 1/6 K coupling for Tb2Ti207.— 



measurements of Ref. [S^- Finally, the diffuse paramag- 
netic neutron scattering patter n^^i^'^i^^i^^ of Tb2Ti2 07 
differs drastically from the experimental spin ice pattern 
(which has been reproduced by Monte Carlo simulations 
of the DSIM^ and its improvements^). This strongly 
suggests that the Ising nature of the locahzed moments 
is not an appropriate description for the magnetism in 
TbaTiaOr, as noted in Ref. H- 

Some important insight into the microscopic nature of 
Tb2Ti207 is provided by a mean field theory for classical 
spins with only a finite Ising anisotropyj^ Specifically, 
Ref. m finds that a toy model in which spins, subject 
to a finite anisotropy and interacting via isotropic ex- 
change and dipole-dipole interactions, captures the main 
features of the experimental paramagnetic diffuse neu- 
tron scattering pattern in Tb2Ti2 07i^ The results of 
Ref. [13 lead one to suspect that the weaker anisotropy 
of the spins in Tb2Ti207, in contrast to those in the 
spin ices, can be attributed to the fact that because the 
ground to first excited crystal field gap is much smaller in 
Tb2Ti2 07, the effect of excited crystal field states can- 
not be ignored. The effects of VCFEs can be studied, 
albeit incompletely, within the random phase approxima- 
tion (RPA). A computation of the RPA diffuse neutron 
scattering intensity in the paramagnetic regime using the 
full crystal field level structure and wavefunctions^ leads 
to results that are in good qualitative agreement with 



experiment,-^ adding weight to the idea that one of the 
effects of VCFEs in Tb2Ti2 07 is to decrease the Ising 
anisotropy of the spins. 

Having identified VCFEs as an important contribution 
to the physics of Tb2Ti207, we look for a way of exam- 
ining the effect of VCFEs on the ground state of per- 
haps the simplest minimal model for Tb2Ti207. An ap- 
proach that is well-suited to this problem is an effective 
Hamiltonian formalism. The low energy theory that is 
obtained within this formalism inhabits a product of two 
dimensional Hilbert spaces — one for each magnetic site 
— spanned by the ground state crystal field doublet. So, 
the effective theory can be written in terms of (pseudo) 
spins one-half. Neglecting VCFEs, the effective Hamilto- 
nian is simply the theory obtained by projecting onto the 
ground state crystal field doublet on each magnetic ion 
which, as we shall sec, is the DSIM of interacting (clas- 
sical) Ising spins i.e. a model in which transverse spin 
fiuctuations are absent The separation of energy scales 
to which we have alluded then allows us to develop a 
perturbation series in the parameter {V) /A — where the 
zeroth order term is the DSIM — and higher order terms 
explicitly incorporate the effect of VCFEs in terms of 
operators acting within the projected Hilbert space. The 
procedure can be written schematically as 



i/(J) =H,f + V 

projection 



perturbation theory 



^^eff (S, 



eff J 



where the bare microscopic Hamiltonian H, depending 
on magnetic moments J through the crystal field iJcf and 
interactions V , is used to derive an effective Hamiltonian 
-ffcff in terms of pseudospins 1/2, Soff. 

One advantage of this approach is that, by decreas- 
ing (y)/A, we can smoothly connect our results to the 
physics of spin ice.— i^*^ A second more practical advan- 
tage is that, since the dimensionality of the relevant 
Hilbert space is reduced, exact diagonalization calcula- 
tions on finite size clusters (albeit small clusters), series 
expansion techniques and the linked cluster method may 
become tractablci^ 

A comparison has previously been made^^ between the 
effective Hamiltonian to lowest order in quantum correc- 
tions, (J'ex)/A, with the crystal field gap A as a free pa- 
rameter and the "high energy" microscopic (bare) model 
from which it was obtained. This involved an exact di- 
agonalization of the two models on a single tetrahedron 
to determine the ground state as a function of A and 
the exchange couplingi^ The result is shown in Fig. [31 
The ground state degeneracies largely coincide over the 
range of parameters explored, which includes the esti- 
mated exchange coupling of Tb2Ti2 07. Most impor- 
tantly, in the singlet region of the phase diagram, the 
ground state of the exact bare microscopic model is a 
nondegenerate superposition of states each satisfying the 
spin ice constraint. In contrast, for the classical dipolar 
ice model with the same exchange coupling, on a sin- 
gle tetrahedron and on a lattice, the ground state is a 
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doubly degenerate all- in/all-out state (see Fig. 14(a)). 



That the full quantum problem favors spin ice-like corre- 
lations at the single tetrahedron level was shown to arise 
from a renormalization of the Ising exchange in the effec- 
tive anisotropic spin- 1/2 Hamiltonian when VCFEs are 
included.— Finally, it was found that the level structure 
from exact diagonalization of the original model on a sin- 
gle tetrahedron is sufhcient to reproduce the main semi- 
quantitative features of the experimental diffuse neutron 
scattering pattern for Tb2Ti207i^ 

The renormalization of the effective nearest neighbor 
Ising exchange by VCFEs such that spin ice correlations 
are energetically preferred over a larger range of the bare 
exchange couplings than would be the case without quan- 
tum corrections shows clearly that quantum effects can 
have a significant effect on the nature of the correlations 
in Tb2Ti207. However, owing to the presence of a long- 
range dipolc-dipole interaction and the fact that VCFEs 
in themselves generate interactions beyond nearest neigh- 
bor, it was not clear on the basis of earlier work^^ whether 
VCFEs would have a significant, or even the same qual- 
itative effect on the Tb2Ti207 correlations when consid- 
ering the full lattice. That is the main problem that we 
resolve in this work. 



B. Scope of the paper 

In this article, we present a more detailed derivation of 
the effective Hamiltonian for Tb2Ti207 than was possible 
in the earlier work^^ owing to lack of space. We also take 
some initial steps beyond the single tetrahedron approx- 
imation by calculating the ground states of the effective 
model assuming that the effective S'eff = 1/2 spins are 
classical spins of fixed length (large S approximation). 
Our main result is shown in Fig. [2] which is discussed 
more fully in Section IV Dl The plot shows the semi- 
classical phase diagram of the effective model on a cubic 
unit cell with periodic boundary conditions as a function 
of the gap A and the isotropic exchange coupling jTcx in 
the microscopic model. When 1/A = 0, all quantum cor- 
rections are suppressed and we recover the limit of the 
dipolar spin ice model (DSIM) with two phases - a state 
with the spin ice rule satisfied on each tetrahedron and 
ordering wavevector 001 (LRSIqoi) and a four-in/four- 
out Ising state (AIAO) for more antiferromagnetic jTox- 
Compared to the dipolar spin ice model ground states, 
the effective model contains one other phase — a q = 
long range ordered spin ice phase (LRSIqoo)- Also, the 
magnetic moments in the LRSIqoo and LRSIooi phases 
are canted away from the local Ising directions as A de- 
creases. The region over which the LRSIooo is the ground 
state forms a wedge, broadening out to lower A until it is 
the only phase found within the explored range of jTex at 
the expense of the antiferromagnetic AIAO phase. There 
are two main physical mechanisms (contributions) to the 
stabilization of the LRSIooo state across the phase di- 
agram. The first is that the effective nearest neighbor 
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FIG. 3: (color online). Figure showing exact diagonaliza- 
tion of a minimal Hamiltonian, H = Hd -I- V, on a single 
tetrahedron. The ground state degeneracy is shown for dif- 
ferent values of the ground-to-first excited crystal field gap A, 
and the bare exchange coupfing for fixed dipolar strength 
© = 0.0315 K relevant to Tb2Ti207. There are two regions: 
one with a singlet ground state, the other with a doubly de- 
generate ground state. The boundary between the two regions 
is marked for the two models considered. For the efi'ective 
Hamiltonian the boundary is marked by circles and for the 
four crystal field state microscopic model on a single tetra- 
hedron described in the main text (based on the crystal field 
Hamiltonian Eq. (|27|) ). the boundary is traced out by squares. 
For the estimated parameters (J'ex,-D,A) for Tb2Ti207, in- 
dicated by a star, the boundaries agree to within ten per- 
cent. The horizontal dashed line shows the phase boundary 
of the classical part (1/A — 0) of the efi'ective Hamiltonian 
between the all-in/all-out doublet configurations and sextet 
(degenerate two-in/two-out) "spin ice" ground states. Within 
the classical description, Tb2Ti207 would be in the doublet 
all-in/all-out state (i.e. above the horizontal dashed line). 
However, when VCFEs are included, the phase boundary is 
shifted towards larger (i.e. more antiferromagnetic) values of 
the bare exchange in such a way that Tb2Ti2 07 "finds itself" 
below the boundary in a singlet ground state. The singlet 
arises because of fluctuations that lift the sixfold degeneracy 
of classical two-in/two-out configurations on a single tetrahe- 
dron. 



Ising coupling becomes more ferromagnetic in character 
as A decreases. However, it does eventually change sign 
as Jcx increases over the entire range of A studied. So the 
second reason for the spreading of a spin ice state across 
the phase diagram as A decreases is due to beyond near- 
est neighbor interactions that arise purely from effective 
VCFEs and which monotonically increase in strength as 
A decreases. 

The outline of the paper is as follows. In Section 
im we introduce some notation and describe the micro- 
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scopic (bare) model for Tb2Ti207 from which the effec- 
tive model is derived. With this in hand, we formulate 
our approach in more detail than in this introduction. 
Section [Till discusses the form and properties of the low- 
est order (classical dipolar spin ice) term in the effective 
Hamiltonian. In Section HVl the quantum corrections to 
this model are enumerated to lowest order in (V) /A and 
we study how the longitudinal (Ising) exchange coupling 
in the dipolar spin ice model (DSIM) is renormalized to 
this order. Having obtained the effective Hamiltonian 
for Tb2Ti207 to lowest order in the 1/A, we treat the 
effective S = 1/2 spins as classical spins and present, in 
SectionlYl the resulting semiclassical ground states. This 
study of the ground states allows us to see how the ef- 
fect of VCFEs is connected to the physics of spin ice and 
also clearly shows that spin ice correlations are present 
even though the bare microscopic exchange coupling jTox 
is antiferromagnetic. 

In other words, geometric frustration in 
the model (Eqs. and 0jj of r&2 Ti2 O7 

emerges from quantum virtual crystal field ex- 
citations (VCFEs) and many-body physics. 

This is the main result of our paper. We discuss 
these results, in Section IVTl in the light of experiments 
on Tb2Ti207 and describe some possible further appli- 
cations of the effective Hamiltonian that we derive for 
Tb2Ti207. Finally, we provide in Appendix 1X1 details of 
the effective Hamiltonian method as a background to the 
main application to Tb2Ti207 described in the remain- 
der of the paper. Appendix [B] contains further details 
behind the calculations presented in Section IIVI and Ap- 
pendix [C] gives some data used to convert between crys- 
tal field parameters for different rare earth pyrochlore 
titanates using a point charge approximation. 

We note here that while our specific focus is on the 
Tb2Ti2 07 pyrochlore magnet, the formalism that we em- 
ploy below could be straightforwardly used to construct 
effective low energy theories for many other frustrated 
rare earth systems where the excited crystal field levels 
have a somewhat larger energy scale than the microscopic 
interactions. 



II. EFFECTIVE HAMILTONIAN 

A. Microscopic (Bare) Model 

The microscopic or bare Hamiltonian for the magnetic 
Tb3+ ions in Tb2Ti2 07 is given by 



H = Hr 



V 



(1) 



where flcf is the crystal field Hamiltonian and V are the 
interactions between the ions. In the remainder of this 
section we explain the form of both terms in some detail. 

The magnetic Tb'^^ ions in Tb2Ti207 are arranged on 
the sites of a pyrochlore lattice. The pyrochlore lattice 



consists of corner-shared tetrahedra which can otherwise 
be thought of as a face-centered cubic (fee) lattice with 
primitive translation vectors for A ~ 1,2,3 and a 
basis of four ions (a = 1, . . . , 4). We follow the same 
labeling of the four sublattice basis vectors as in Ref. I&j . 
It is useful to introduce a coordinate system on each of 
the four sublattices with local z° unit vector along the 
local cubic [111] direction. The sublattice basis vectors 
and local Cartesian x"^, y" and z° directions are given 
in Table IH Below, we also make use of rotation matri- 
ces w° ^ (the elements of which are contained in Table |T| 
which achieve a passive transformation that takes the lo- 
cal sublattice coordinate system for sublattice a into the 
global Cartesian laboratory axes. 

Spin-orbit coupling within the relevant localized 4/ 
levels of the Tb''^ ions leaves total angular momentum J 
as a good quantum number with J = 6. The local envi- 
ronment about each Tb'^^ ion is responsible for breaking 
the 2J -I- 1 degeneracy. Its effect can be computed from 
a crystal field Hamiltonian, Hcf, which is constrained by 
symmetry to take the form^i^i^ 

TJcf = ^2 0^(^ a) + BlOlii, a) + 0|(i, a) 

i,a 

+ i?°06"(*, a) + BlOUl, a) + BtOl{i, a). (2) 

The magnetic ions are labeled by an fee site i and a sub- 
lattice index a. Expressions for the operators O™ in 
terms of the local angular momentum components can 
be found, for example, in HutchingSi^ The crystal field 
in Tb2Ti207 has been studied in Refs. [10 andllsl result- 
ing in somewhat differing estimates for the parameters 
S™. In the following, all quantitative results that we 
present for Tb2Ti2 07 were obtained using crystal field 
parameters for Ho2Ti207. obtained from inelastic neu- 
tron scattering in Ref. [43, which have been rescaled to 
the Tb2Ti207 parameters according to 



Tb 



({Si) 



Tb 



Ho 



Tb 



Ho 



Ho- 



(3) 



Here, the Si are Stevens factors.— These and the radial 
expectation values (r™) for the rare earth ions^ can be 
found in Appendix [C] We have checked that using the 
crystal field parameters of Ref. \^ instead leads to results 
that are in fairly close quantitative agreement with those 
obtained using the rescaled parameters from Eq. 

The crystal field Hamiltonian, i?cf, can be diagonal- 
ized numerically exactly; the eigenvalues are En and the 
eigenstates |n) for n = 1, . . . , 13, which we implicitly ar- 
range in order of increasing energy. One finds a level 
structure that includes a ground state and a first excited 
state that are both doubly degenerate . The splitting, 
A, between the ground and first excited states is about 
18.6 K)^i^ which is much smaller than the correspond- 
ing gap in the spin ices (for example, the gap in Ho2Ti207 
is about 230 K— ). It is the smallness of this value of A 
compared to V for Tb2Ti207 and the possibility of ad- 
mixing between the ground state and excited state crystal 
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TABLE I: Basis of four magnetic ions on a pyrochlore indexed by position vectors r". The local [111] direction on each sublattice 
is z". The edge length of the cubic unit cell is a. The rotation matrix u^^ takes the form (x", y", z")^^ in which the vector 
components are placed in the matrix columns. In the main text, we make use of vectors = (1, 0, 0),n'' = (0, 1, 0),n' = (0, 0, 1) 
in the laboratory coordinate system. 



Sublattice 



1 (a/4)(0,0,0) (l/^)(-l,-l,2) (1/V2)(1,-1,0) (1/V3)(1, 1, 1) 

2 (a/4)(l,l,0) (1/V6)(1,1,2) (1/V2)(-1, 1, 0) (l/^/3)(-l, -1, 1) 

3 (a/4)(l,0,l) (l/^)(l,-l,-2) (i/y2)(-l,-l,0) (1/V3)(-1, 1, -1) 

4 (a/4)(0,l,l) (l/^)(-l,l,-2) (l/x/2)(l,l,0) (1/V3)(1, -1, -1) 




18.6 K 



FIG. 4: (color online) . Figure indicating the four lowest levels 
of the crystal field spectrum (not to scale). The splitting 
between the ground state doublet and the first excited state 
is called A. The ground state and the first excited state are 
doublets.— The two other excited states are singlets.— 



field levels that are at the root of all the phenomenology 
that we explore in the rest of this paper. Fig. |4] shows 
the level structure of the crystal field spectrum for the 
four lowest levels determined on the basis of an exact 
diagonalization of Eq. ([2]) . 

We emphasize two features of this spectrum that will 
be important later on. First of all, let us write down the 
time reversal properties of the eigenstates, |n). Let |n) 
be written as a linear combination of the eigenstates of 
J, denoted |J,M), 



|n) = 5]cf|J,M). 



M 

Time reversal invariance requires that the coefficients are 
related to one another by c*^ = (— )''~*^c~*^.— Secondly, 
it is possible to interpret the non-interacting single ion 
angular momenta as Ising-like at low energies, as was 
done in Ref. [1^. This is because, at sufficiently low ener- 
gies, thermal occupation of excited crystal field levels is 
negligible and one can focus on the ground state doublet. 
The ground state doublet states, |1) and |2), have 

(1|P|1)=-(2|P|2)^(P) (4) 
as the only nonvanishing matrix elements, where the tilde 



indicates that the z axis is taken along the local [111] di- 
rection appropriate to each magnetic ion (see Table 
So, this doublet considered on its own has nonzero an- 
gular momentum expectation values only along one axis 
with vanishing transition matrix elements (1|J^|2) = 0. 

The interactions between the angular momenta, V = 
Hex + Hdd, are taken to be nearest neighbor isotropic 
exchange H^x and dipole-dipole interactions, Hdd- 

Hex *!Tcx ^ ^ Ji,a ' 



E 



-3 



The notation R?, is short for R? 



-Rj withR? 



(5) 



R, 



and = 3.59A = a\/2/4 (where a is the edge length of 
the conventional cubic unit cell) is the distance between 
neighboring magnetic ionsi^ Here, we employ the con- 
vention that i/ex > is antiferromagnetic and J^x < is 
ferromagnetic. This is the simplest Hamiltonian consis- 
tent with the nonvanishing Tb''^ dipole-dipole coupling, 
V = {^io/'iTT)igjHB f/rl^, = 0.0315 K with the Lande 
factor, gj = 3/2 and with the negative Curie- Weiss tem- 
perature ^cw = — 14 Ki^ The exchange coupling jTox has 
been estimated from 9cw for Tb2Ti2 07 and 6cw for the 
diluted compound (Yo.98Tbo.o2)2Ti207 — to be about 
0.17 K, while a fit in ReL [H gives a value for jTcx = 0.083 
K that is significantly less antiferromagnetici^ 

In summary, our bare microscopic model for Tb2Ti207 
consists of three terms: the crystal field Hamiltonian 
i?cf, an isotropic exchange H^x with an antiferromag- 
netic coupling and a dipole-dipole interaction, H^d'^ 
An extension of the present work could include (i) 
bare exchange couplings beyond nearest neighbors, (ii) 
anisotropic nearest neighbor exchange as described in 
Ref. m and (iii) direct or virtual (phonon- mediated) mul- 
tipolar interactionsi^ 



B. Route to an effective Hamiltonian 

If we were able to ignore the excited crystal field lev- 
els in Tb2Ti207, the angular momenta could be treated 
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as classical Ising spins^ii^ because the only nonvanish- 
ing matrix elements of the angular momentum are those 
in Eq. ([4]).— However, for reasons outlined in the Intro- 
duction, this is not a good approximation for this ma- 
terial. The interactions between the angular momenta 
induce VCFEs that admix excited crystal field wavefunc- 
tions into the space spanned by the non-interacting crys- 
tal field doublets with the consequence that the mag- 
netic moments behave much less anisotropically than 
one would expect on the basis of the [111] Ising- like 
ground state crystal field doublet. These quantum fluc- 
tuations can be treated perturbatively because there is 
a small dimensionless parameter (l/)/A, where {V) ^ 
0(max(j7ex, 2^))- To lowest order in such a perturbation 
theory, and in a low energy effective model, the spins 
should be perfectly Ising-likc and hence we recover the 
DSIM. We now proceed to make these ideas more con- 
crete. 

Because we seek a Hamiltonian operating within a low 
energy subspace, we need a projection operator onto the 
non-interacting single ion crystal field ground states. For 
a single ion at the site specified by indices i,a, the pro- 
jection is accomplished by 

r{i,a)^\h^a){l^.a\ + \2^,a){2^^a\. 

This operator satisfies the conditions V'^{i,a) = V{i,a) 
and Hermiticity. With moments on all the sites of the 
lattice, the projector is T' = Jlia^C*''^)- subspace 
of the full Hilbert space selected by the projector will be 
called the model space, 9Jl = Yl^(^i 2ti,a, from now on. 
The Hilbert space Tli^a is defined as the space spanned 
by states |li,o) and \2i^a) on site (i,a). 
The spin-spin interaction 



V = H, 



H, 



dd 



(6) 



is to be treated as a perturbation. Because the per- 
turbation V is "small" compared to the difference be- 
tween the ground and first excited crystal field energies 
A, Hci = Hq, we expect that on a crystal of N sites, the 
2^ lowest energy eigenstates of H lie mainly within SDt 
because the admixing of excited crystal field wavefunc- 
tions into the model space is a small effect. Our effective 
Hamiltonian will be defined in such a way that its eigen- 
states live entirely within 971 while its eigenvalues exactly 
correspond to the 2^ lowest energy eigenvalues of the ex- 
act Hamiltonian, H. The 2^ lowest energy eigenstates 
of H mainly lie within 971 in the sense that the rotation 
of exact states out of the model space is determined by 
the relatively small perturbation {V) / A. 

In practice, the exact eigenvalues can be approximated 
by carrying out perturbation theory in the construction 
of the effective Hamiltonian i?off • After some work, that 
is briefly laid out in Appendix \^ one finds that the ef- 
fective Hamiltonian can be written as&^ 



The operator TZ — the resolvent operator— is given by 



En — Ef 



(8) 



VHoV + rvv + vvnvv 



(7) 



where, for a finite crystal of N sites. Eg is N times the 
energy of the degenerate ground state crystal field lev- 
els Eq. The numerator of each term in the resolvent is 
a projector onto a space orthogonal to 971 — a product 
of crystal field operators \P){P\ = Yl^ where the 

product is taken over all sites of the lattice with at least 
one such operator having n > 2 (i.e. belonging to the 
group of excited crystal field states); this is the meaning 
of the notation |P) ^ 971 in the summation index of Eq. 
(|S]) . The third term on the right-hand-side of Eq. ([7]) is 
the lowest order term in the perturbation series to include 
the effects of crystal field states outside the model space. 
This term is therefore the lowest order contribution of 
the VCFEs that we have referred to above. 

Equation ([7]) makes no reference to a particular model. 
In Sections IIIII and IIVI we develop the terms in the 
effective Hamiltonian for the model H = Hq + V = 
Hcf + H,id + Hex of Tb2Ti2 07 described in Section HT^ 
Section IIIII is devoted to the lowest order, or classi- 
cal, term VHV. Section HVl enumerates the lowest or- 
der terms generated by VCFEs, relating each underly- 
ing class of terms that originate from VHTZHV to spe- 
cific virtual excitation channels. Higher order corrections 
than VHTZHV are computationally difficult to determine 
mainly because of the presence of the long-range dipole 
interactions iJdd- See Ref.[6lfor a model on a pyrochlore 
for which degenerate perturbation theory can be carried 
out to much higher order than is done is this work. 

To spare readers the details of this rather technical 
derivation if they so choose, we include a short summary 
(Section IIVF[) of the form of the low energy model for 
Tb2Ti207. Finally, in Section [IV Gl we summarize some 
results that have been obtained from the effective Hamil- 
tonian which have already appeared in the literaturei^i^ 
All in all, we shall see that the DSIM couplings arc renor- 
malized by VCFEs and that effective anisotropic spin- 
spin couplings appear in addition to the Ising interactions 
of the DSIM. In other words, the effective theory allows 
for fluctuations of the moments perpendicular to the lo- 
cal z° axes. We shall study the variation of the effective 
couplings in Hcs as jTcx is varied. This information will 
be useful in the interpretation of the semiclassical ground 
states of the effective model (Section |V| and hence in as- 
sessing the effects of VCFEs on the physics of Tb2Ti207. 



III. CLASSICAL PART OF H^s 

A. [Ill] Ising model for Tb2Ti207 

In this subsection, we consider the (lowest order) term 
VHVinEq. ([7]). The effective Hamiltonian derived from 
H for Tb2Ti207 can be rendered in the form of a spin 
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one-half model by rewriting the model space operators 
in Eq. ([7]) in terms of Pauli matrices. This is possible 
because the model space, in our case, is a direct product 
of two dimensional Hilbert spaces spanned by the ground 
state crystal field doublet. The correspondence between 
Pauli matrices and operators on the crystal field ground 
state: 



a- = |l)(2| + |2)(l| 
a^ = -*(|l)(2|-|2)(l|) 
a^ = |l)(l|-|2)(2| 



(9) 
(10) 
(11) 



together with the unit operator I = |1)(1| + |2)(2|. Note, 
however, that despite the fact they do satisfy the com- 
mutation rules 



^1 



where Ca^-y is the Levi-Civita symbol, the ct" do not swap 
sign under time reversal so they are not true angular mo- 
mentum operators. For this reason, we shall call them 
pseudospins or effective spins. It is helpful for later sec- 
tions to give their properties under T, the time reversal 
transformation: 



T : CT?' ^ CT?- 



(12) 
(13) 
(14) 



because 1 : |1(2)) ^ |2(1)). 

If we apply the projector to the full Hamiltonian to 
obtain VHV, we find that the crystal field part Hd be- 
comes EoJ2i a^i.a 'with El = E2 = Eq. From now on, 
we omit this constant energy shift. To project the in- 
teraction part VW, we write the angular momentum 
components in the local coordinate system with local z" 
axes in the directions given in Table [J J"^ = w^^J^^^. 
All operator components that refer to the local coordi- 
nate systems are labeled with a tilde. Also, when it is 
not important to distinguish different sublattices, we ab- 
breviate (i, a) with the site index /. We add further nu- 
merical subscripts to / to label different sites. With this 
notation, the projector acting on J| gives 



(^l)(|l/)(l/| 



\2j){2j\) = {J})af 



where (J|) = (IjJ^l). Owing to (1|J^|2) = 0, all matrix 
elements of the other angular momentum components 
vanish. So, the isotropic exchange H^x becomes 



i^lassical ^ ^ (z 



and the dipolc-dipole interaction becomes 



(15) 



1 



{i,a;j,b) 



(z'^.z") Jz°-R^'')(z'' 



-3- 



\.b- 



(16) 



The renormalized, or effective, exchange and dipolc- 
dipole couphngs are^ respectively, Jciassicai = Jcx(./^)^ 
and ^classical = V{.F-)\ HusM - ViH,x + Hm)V is 
the celebrated DSIM.i>^>^>iiii^>^>^ It is a classical (local 
Ising) model because all the terms mutually commute as 
they solely consist of trf ^ operators. This model exhibits 
two different ground states depending on the ratio of the 
exchange to the dipolar coupling; these are shown in the 
inset of Fig. [H When Je^/V > 4.525, the ground state 
has ordering wavevcctor q = with the spins on a single 
tetrahedron in the |1) state or the |2) state — the all- 
in/all-out phase4^ When Jc^/V < 4.525, the ordering 
wavevector of the ground state is (0, 0,27r/a) and each 
tetrahedron has spins satisfying the two- in/two-out ice 
rule; we refer to this state as the LRSIqoi phase,— i^i^ 
with one of the domains shown in Fig. [l}^ Above a 
nonzero critical temperature, the LRSIqoi phase gives 
wayi to a spin ice state with no conventional long-range 
order (Fig.©. 

Formally speaking, the spin ice state is a collective 
paramagnetic stated — a classical spin liquid of sorts. 
That the DSIM has proved to be a good model for spin 
ice materials is largely due to the substantial gap A be- 
tween the crystal field ground state doublet and first ex- 
cited state which results in a roughly 1/A suppression 
of VCFEsi^ This model is not a good description for 
Tb2Ti207. Indeed, if we consider the estimated cou- 
plings given in Section III Al we find Je^jT) ~ 5.4 (re- 
calling J7cx ~ 0.17 K and T) k, 0.0315 K as stated in Sec- 
tioning, which would put Tb2Ti207 in the all-in/all- 
out phase with a critical temperature into this phase 
from the paramagnetic phase at ^ 0.5 K (see vertical 
dashed line in the inset to Fig. [2]) 4^ This is in contra- 
diction with neutron scattering experiments which find 
no magnetic Bragg peaks in zero field 1^ If we allow for 
inaccuracies in the estimate of jToxr^'^ such that a clas- 
sical, dipolar spin ice state is implied by the coupling, we 
find various properties of spin ices that are not compat- 
ible with those of Tb2Ti207. Some of these confiicting 
properties — the diffuse neutron scattering pattern and 
differing spin anisotropics — were discussed in the Intro- 
duction. Therefore, in the next section, we investigate 
what happens when A is small enough that the lowest 
order fluctuation term VHTZHV in Eq. ([7]) becomes im- 
portant. 



B. Exchange convention 

In Eq. ([5]), we use the opposite sign convention for the 
exchange coupling to the one used in Refs. |4l|, and 
l48li^ The convention in these works is to include a minus 
sign in front of the exchange coupling in contrast to our 
Eq. ([5]). In this article, in the global coordinate system, 
antiferromagnetic corresponds to jTex > 0. 

A warning must be made regarding the convention 
within the local coordinate system. In rotating to the 
local system, geometrical factors appear in front of the 
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couplings. For example, as shown in Section ITTll the lo- 
cal Ising exchange part of the coupling j7exJi,a • Jj.fc is 
Jo-^iif^ ■ z^)Ji,a ■ Jjb where (z° • z^) = —1/3 arises from 
the fact that the local z axes are not collinear. In the 
following pages, we adopt the simplifying scheme of ab- 
sorbing the geometrical factors into the couplings. In 
doing so, it will be useful to describe how to go from 
the sign of the local effective Ising coupling, J^^ , in front 
of J^^a^a^ to the type of order that is energetically fa- 
vored by the coupling. Thus, when the local coupling is 
said to be ferromagnetic, J^^ is negative and the Ising 
components of the spins prefer to lie in an all-in/all-out 
configuration. When, instead, J^^ is positive, it is said 
to be antiferromagnetic and the local Ising components 
are frustrated, leading to a spin ice configuration on each 
tetrahedron. 



IV. LOWEST ORDER QUANTUM 
FLUCTUATIONS 

A. General Considerations 

In this section we present a derivation of the quantum 
terms VVTZVV in the effective Hamiltonian. We refer 
those readers interested only in the results of this techni- 
cal derivation to Sect ion HV Fl We begin by introducing a 
little more notation to describe the structure of the term 
VVTZVV which we shall refer to as H^g . We write the 
interaction term V in the form 

I1J2 a,P 

= Y.T.^tiAJl (17) 

where, in the second line, we have absorbed the rotation 
matrices into the definition of JC. When the spins 
interact via nearest neighbor isotropic exchange and long- 
rangcd dipole-dipole interactions as in Eq. ([5]), we have 
for /C: 

with unit vectors n" for a ~ x,y,z m the laboratory 
x,y,z directions respectively (see Tabic |l|. The prefac- 
tors of one-half cure the double counting of pairs in Eq. 

The model space M basis states are products of ground 
state doublet states |1/) and \2j) over all lattice sites / 
while excited crystal field states are denoted \Wi) for 
ly = 3, . . . , 13 on each site /. The state |P) in Eq. (O is 
a direct product of crystal field states on different sites 
with the condition that at least one of the states in \P) 
lies outside the ground state crystal field doublet; in other 



words, at least one Tb^+ ion must be virtually excited in 
a state \n) with n > 3. With this notation in hand, we 

(2) 

write the quantum term H^^f' = VVTZVV as 

E E T'{'^tiJTA)'^{nfijiJi)'p. (19) 

There are a few observations that we can make from 
Eq. that identify classes of nonvanishing terms. 

Suppose we choose magnetic sites Ip on the pyrochlore 
lattice for p = 1,2,3,4 in Ea. p9)) . Then, when we eval- 
uate Eg. p^ for all other sites, we obtain unit operators 
|l/™>(l/™l + |2/™)(2/„| for all sites with m 7^ 1,2,3,4. 
This follows because the resolvent operator TZ and pro- 
jection operators V are diagonal on each site. In the 
following, we do not write out all these unit operators ex- 
plicitly. A second observation is that when we consider 
a term with magnetic sites Ip [p = 1,2,3,4) all differ- 
ent, we find that such a term vanishes. The reason for 
this is that the resolvent and angular momentum opera- 
tors are sandwiched by projectors into the model space. 
That way, a virtual excitation induced, for example, by 
J/3 in the bilinear operator must be "de-excited" 

by an angular momentum operator in the other bilinear 
operator /C/j/j, (with Ii = /a, for example). If all Ip are 
different there can be no virtual excitations and, because 
the resolvent operator is orthogonal to the model space 
states, such terms must vanish. 

Having found those terms that must always vanish, we 
now divide all the potentially nonvanishing terms into 
three classes that we shall analyze in turn in the next 
three subsections. 

CASE A The first class of terms has two groups (/i,/2) 
and (/a, I4) of sites, with exactly one site in the 
first group in common with a site in the second 
group. In this case, the ion on the common site 
must be virtually excited and de-excited, and 
the other two ions remain in their ground dou- 
blets. This is because the resolvent operator 
demands that there be some virtual excitations 
and that the projectors require that any virtual 
excitation must be dc-excitcd. So. only when 
two angular momentum operators (one in each 
V operator of VVTZVV) belong to a given site 
can that site be virtually excited. 

CASE B This class of terms has identical pairs (/i,/2) 
and {13,14) regardless of label ordering, but 
with only one single ion (Ji or I2) that is vir- 
tually excited. 

CASE C Finally, we shall consider the case where (/i , I2) 
and (I^jli) are identical pairs and where both 
ions are virtually excited. 

The virtual excitations belonging to each of these three 
cases are illustrated in Fig. O 

It will be convenient, while considering the possibilities 
enumerated above, to make use of the following explicit 



10 



/ \ 



part and VH^dR-UddP as the dipolc-dipolc part. 



\ 

(a) An example of Case A with I2 = I4 and Ii Is- 
Only the ion on site I2 is virtually excited. The 
two pairs of sites Ii , I2 and I2 , I3 are not 
restricted to be nearest neighbors because the 
bare Hamiltonian has long-ranged dipole-dipole 
interactions. 



/ \ 



I2=l4 

I4 with virtual 



I,= l3 

(b) Case B with /i = J3 and I2 

excitations only on site 12- There is no virtual 
excitation on site Ii. I\ and I2 need not be n- 
earcst neighbors because K includes an intera- 
ction between dipole moments. 



/ \ / \ 



I, = I, I2=l4 

(c) Case C with /i = Ii and I2 = Ia- Ions on both 
sites arc virtually excited. Similarly to Cases A 
and B, Ii and I2 need not be nearest neighbors. 

FIG. 5: (color online). Figure illustrating the virtual exci- 
tations distinguishing three classes of terms in the effective 
Hamiltonian which are enumerated and described in the main 
text. The arrows show virtual excitations and de-excitations 
within the lowest-lying pair of crystal field doublets belonging 
to the ion on the labeled site. Sites with a blue circle over 
the ground state doublet indicate that the ion on that site 
remains in its original state within the ground state crystal 
field doublet. 



(2) 

decomposition of the quantum term H^^ : 

vvnvr = 7'i^ex7^i^exP + {VH^^nH^aV 

+VHMnH,^V) + VH^nH^dV. (20) 

We refer to VHcyJZHc^V as the exchange-exchange part, 
{VHcxTZHm'P + 'PHmT^Hcx'P) as the exchange-dipolc 



B. 



Case A 



We will show that the situation in Case A described 
above leads to (i) effective Hamiltonian bilinear interac- 
tions between the local z components of the spins and 
also to (ii) three-body interactions of the form o-f^o-f^o-f^ 
where a = x or y, but not z. 

We write the bilinear operator on the left-hand-side 
of Eq. as and the other bilinear as 

'^hiJf2Jh ^ith h ^ h ih = h, see Fig. with 
all angular momentum components referred to the local 
coordinate system. As we discussed above, the contribu- 
tion of all the other sites gives identity operators for each 
site. Omitting these unit operators, we are left with 



E EE^(-i'™2,m3)i^;5^j,v^: 

aJ3,p,a nip W 

1^4,71 , , 7713 j^) {1114 , Wi^ , TO3 J3 I 

X 

Eq — Ew 

X KP^.JfJinm^, m,,me) (21) 

with 

7^(7711,7712,7713) = |r77lJ^,7772,72,"73,73)("^l,7i,7772,72,nT-3,73l 
7^(7714,7775,777.6) = I7774. 7^,7775,72, 7776,73) (7774 Jj, 7775, 72,7776,73 I . 

(22) 

Ew is the energy of an excited crystal field state on a 
single ion. Here, the angular momentum components 
are expressed in their respective local coordinate systems 
with local z axes given in Table |T1 the rotation matrices 
having been absorbed implicitly into ^7^72 ■ The integers 
rUp run over 1 and 2 with the states lying within DJlj on 
their respective sites. We factor out the part for site Ii: 
Smi,m4 l'^i)("^i|"^/il™4)(7774|. Recalling the property, 

Eq. ^ and the mapping in Eq. (fTTj) . we obtain { J^)aj^ . 
We reach the same result for the sum over states 7773 and 
7776 on site I3. Equation (|2ip then simplifies to 



TP 



E El"^2)(7775| 



(7772|J,||l^)(M^|j;j7775) 

Ea~Ew 



(23) 



in the local coordinate system where we have dropped the 
I2 site labels from the state vectors enclosed by brackets. 
After summing over the excited states \ W), and render- 
ing the sum of operators in terms of Pauli matrices, we 
find an Ising interaction term cr|^ af^ and three-body op- 
erators (t|^(Tj^(t|^ with a = x,y and a ^ z since, if a 
were to equal z, the three-body term would violate time 
reversal invariance. The sum over virtual excited state 



and the subsequent rendering in terms of Pauli operators 
is discussed in some more detail in Appendix [Bl 

We have reduced the most general three ion terms in 
H^^ (case A) to interactions between pseudospins one- 
half but we have not made any assumptions yet about the 
form of the interactions K^'jfi^- In the following, we shall 
consider the four terms of Eq. ([20|) in turn within Case 
A. These terms determine the spatial range of the result- 
ing effective interactions between the pseudospins and 
are obtained by distinguishing the exchange and dipolar 
parts of K in Eq. ((23l) . 

Exchange-exchange part The exchange-exchange part 
(referring to the first term of Eq. (|20|) ) — which is nothing 
more than Eq. for K with I? = — is a short-range, 
but not strictly nearest neighbor, effective interaction. 
If Ii, I2 and I3 all lie on the same tetrahedron in the 
lattice, then the Ising interaction acts between nearest 
neighbors and, for given Ii and I3, there are two choices 
for the position of the "mediating ion" I2 as shown in 
Fig. [HI The thick lines in this figure join the Ii and 
I3 ions via two different choices for the ion I2 and gen- 
erate, together with three-body interactions connecting 
/i, I2 and I3, a nearest neighbor effective Ising interac- 
tion between Ii and I3. As we describe in more detail 
later on, this Ising interaction renormalizes the i/dassicai 
exchange, defined in Eq. (|15p . Within a single tetrahe- 
dron, the three-body interactions couple all three pseu- 
dospins along each of the paths in Fig. [B] There are two 
other exchange-exchange pseudospin terms arising from 
Eq. — those for which the interaction extends fur- 
ther than a single tetrahedron. The ions at the endpoints 
of the lines on the left-hand-side and right-hand-side of 
the Fig. [7]are coupled, respectively, by effective Ising-hke 
second nearest neighbor interactions and effective third 
nearest neighbor Ising exchange. The ions along each 
line - those at the endpoints and the one at the center 
of each line - interact also via effective three-body inter- 
actions. The dependence of the effective couplings for 
second and third nearest neighbor Ising interactions, on 
the bare exchange coupling J7ex is shown in Fig. [8] (the 
dipole-dipole couphng V being set equal to zero). The 
second nearest neighbor effective coupling is antiferro- 
magnetic, and the third nearest neighbor effective cou- 
pling is ferromagnetic when the bare exchange coupling 
is antiferromagnetic (J7cx > 0). See Section fill B I for the 
convention on the exchange that we use in this paper. 
We note that there are two distinct types of third near- 
est neighbors on the pyrochlore lattice (see, for example 
Fig. 1 of Ref. Ill or Fig. 2 in Ref. I?!) and that only 

one type — those connected by one mediating ion, or two 

(2) 

edges, along the lattice — appear in H^ff . Third nearest 
neighbours of the other type have sites that are connected 
via three edges along the lattice (as shown by the dot- 
ted line of Fig. [7|) and hence couplings between ions on 
these sites would require two mediating ions to appear 
in the effective Hamiltonian. But, to the (second) order 
of perturbation theory that we are considering, there is a 
maximum of one mediating ion so such couplings do not 
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FIG. 6: (color online). In deriving H^J^^ = VVTZVP, where V 
is a sum of bilinear interactions, we consider a single pairwise 
interaction in the right-hand interaction V on sites I2 and 
and a pairwise interaction in the left-hand V between 7i and 
l2. This choice of terms in -ff^g is referred to as Case A in 
the main text. Because Ji 7^/3, and because the operators 
on these sites are sandwiched between P projectors, the only 
virtually excited site is I2 while the other two sites remain in 
their (noninteracting) crystal field ground state. As we show 
in Appendix [B] one obtains effective o"^ operators on site Ji 
and I3. There are two possibilities for the effective operator 
on site I2 after calculation: it could be a unit operator leaving 
an Ising coupling between sites Ji and I3, or it could give rise 
to a transverse operator af^ or a^^, generating an effective 
three-body term connecting sites Ii, I2 and I3. The figure 
shows a single tetrahedron in a pyrochlore lattice with the 
two ways of joining sites Ii and I3. The total effective Ising 
exchange for this Case A between ions Ii and J3, Jf^i^ is the 
sum of the Ising exchange terms corresponding to each path 
in the figure. 



appear m H^g . 

Having discussed the exchange-exchange part, we 
switch on the dipolar interaction. In doing so, the in- 
teractions become anisotropic even in the bare Hamilto- 
nian. Nevertheless, in Case A, the types of couplings that 
arise in the presence of the dipolar coupling are the same 
as those arising in the exchange-exchange case, the only 
difference being in the range over which the couplings 
act. 

Exchange-dipole part The exchange-dipole part of the 
Class A interactions consists of those terms in Eq. ([25)1 
with I?ox = between, say Ii and I2 and ^/ox = be- 
tween I2 and I3. There are then long-range Ising-like 
interactions between Ji and I3 although the exchange 
Hamiltonian constrains two sites, Ii and /2, in this case, 
to be nearest neighbors. The bare microscopic dipole 
interaction acts between sites I2 and I3 and decays as 
jR/j/gl"'^. Overall the Case B effective interactions de- 
cay as |R/j/3 1""^ at long distances just as the bare dipole 
interactions do on their own. The same is true, by 
symmetry, of the Class A interactions belonging to the 
dipole-exchange term in Eq. (PO)) . There are also three- 
body interactions originating from the exchange-dipole 
and dipole-exchange parts of Eq. ([^0]) where two of the 
spins operators must lie on nearest neighbor sites. 

Dipole-dipole part Finally, the terms in Class A be- 
longing to the dipole-dipole term of Eq. (|^D|) (jTcx — 
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FIG. 7: (color online). Part of a chain of tetrahedra in a 
pyrochlore lattice. For case A (with only exchange interac- 
tions in the bare Hamiltonian) , the quantum part, H^:^ — 
VVTV/V, of the effective Hamiltonian has nonvanishing cou- 
plings for interactions between the sites connected by the 
thick lines. The (green) curve on the left-hand-side indicates 
that Ising exchange acts between second nearest neighbors. 
The (blue) curve on the right-hand-side indicates that those 
third nearest neighbors lying along chains of sites are cou- 
pled by Ising exchange. The red dashed curve indicates a 
second type of third nearest neighbor for which no effective 
coupling is generated as explained in the main text. Three- 
body interactions also arise that couple all three ions along 
each of the paths joined in the figure. These further neighbor 
couplings receive a further contribution when the dipolar cou- 
pling is switched on. Also, in the presence of the dipolar inter- 
action, there are off-diagonal effective couplings of the form 
Jj^i^^j-^S'j^ with fi ^ z, v ^ z in addition to an Ising coupling 
that renormalizes the classical (Ising, first order) DSIM VVP 
term. 



in Eq. (US])) have a range that is a product of two 
dipole interactions so the overall two-body interaction 
between sites Ii and ly, decays as a function of I2 as 
|I^-fi/2l~^|R--r2-f3l~'^ before summing over. 

In summary, we have learned in this section, in the dis- 
cussion following Eq. and in Appendix [Bl that the 
effective pseudospin operator on the "connecting site" I2 
involved in the virtual excitation process can be a unit 
operator so that the resulting nontrivial effective interac- 
tion is a bilinear of Ising pseudospin operators between Ji 
and /a. In this circumstance, there are contributions to 
this interaction for I2 positions at arbitrarily large dis- 
tances from Ii and /a. Together with two body Ising 
interactions, there are also three spin long-range interac- 
tions in this group of terms following from the argument 
given above. There are no constraints on the positions of 
the coupled sites Ji, I2 and /a in the lattice. 



C. Case B 

Now, referring to Eq. we impose the constraint 

/3 = /i and Ia = I2 and, without loss of generality, sup- 
pose that only the ion on site I2 is virtually excited (see 
Fig. 5(b) ). After computing Eq. with the aforemen- 
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FIG. 8: (color online). Further neighbor Ising interactions 
generated by quantum fluctuations when = 0. The upper 
curve is the second nearest neighbor Ising coupling (between 
spins with sublattice labels 1 and 3). It is antiferromagnetic 
in nature. The lower curve is the ferromagnetic third nearest 
neighbor interaction (between two spins with sublattice label 
1 — see straight right-hand blue line in Fig. [7]). When 7? 7^ 0, 
these couplings are renormalized by other contributions from 
Class A and Class C terms. Since couplings are referred to 
operators in the local coordinate system, they contain geo- 
metrical factors from the rotation matrices (see Sections III Al 
andEHl). 



constant shift in energy. As a first step, we write out 
Eq. ([TT]) with the unit operators on all sites but /i and 
I2 factored out: 



'/3 



a,/3,p,cr nup W 



Eq — Ew 



KZ,Jf,J?,Vimi,m,), (24) 



tioncd constraint, we obtain a result that is proportional 
to the unit operator. That is, this case gives rise to a 



where the notation T' (mi, 7712) stands for operator 
I"^i,/i,™2,72)('^i,/i7™2,72l- Because only the ion on site 
I2 is virtually excited, the operator Jf^ acts entirely on 
states within 971/^ so, referring to the matrix elements in 
Eq. ^ we obtain, for ion /i, 

|1/.)(1/J(1/J^IJ1/.)' + |2/,)(2/J(2,JJ|j2,j2 

^{J},)%, (25) 

— the unit operator acting on /i, which we can omit in 
the following. When the sum is performed over excited 
crystal field states on ion I2 , the resulting operators map 
to the unit operator I/^ and Pauli matrices af^ and a^^ 

— all the resulting operators being consistent with time 
reversal. This calculation is performed along the lines 
described in Appendix [Bl 



13 



It would thus seem that, together with the constant en- 
ergy shift, there are effective nontrivial single-site trans- 
verse field operators in the effective Hamiltonian of the 
form and . However, these effective transverse field 
terms cancel on any given lattice site when one sums 
the contributions to these single-site operators coming 
from the bare microscopic pairwise interactions in V . 
Without going into the details of the sum over differ- 
ent contributions to site J2 we see that there must be 
such a cancellation because neither the original model nor 
the effective Hamiltonian formalism distinguishes partic- 
ular directions (as opposed to particular axes) on in- 
dividual sites. When the Tb^+ ions are randomly di- 
luted with nonmagnetic ions, as in (TbpYi_p)2Ti207,— 
there is no longer perfect cancellation of the effective 
single site operators. These effective fields have the ef- 
fect of splitting the degenerate |1), |2) doublet on each 
ion for which the cancellation does not occur. So, in 
fact, the low energy effective theory of the diluted com- 
pound (TbpYi_p)2Ti207 would be somewhat different 
than that of the pure Tb2Ti207 by admitting effective 
random transverse fields. The possible generation of ef- 
fective random transverse fields generated by dilution 
in (TbpYi_p)2Ti207 had previously been proposed in 
Ref. Its]. In the remainder of the article, we shall assume 
that the magnetic ions are not diluted. 



D. Case C 



The class of terms where the two ions I\ and I2 are 



the cutoff coming from the bare exchange ensures van- 
ishing of the effective interaction beyond nearest neigh- 
bor for Case C. This means that the only interactions in 
Case C extending beyond nearest neighbors come from 
the dipole-dipole part of Eq. ((20|) and, because there are 



both virtually excited (see Fig. 5(c) I is the most com- 
plicated of the three cases A,B and C in the sense that 
all two body terms consistent with time reversal and the 
lattice symmetries can and do arise. Because the dipolar 
coupling is nonzero, long range effective interactions ap- 
pear in -ffeff- The calculation of the types of terms and 
their couplings in Case C is most easily accomplished 
by the projection method given in Appendix |B] In this 
section then, we give only the results of this calculation 
— the means of calculation having been outlined in the 
discussion of Section IIVBI and in Appendix [B) As with 
the terms in Case B, the net single-site "fields", (t^ and 
cancel, leaving the Ising interaction , and the 

transverse exchange interactions af^d^^ where each of a 
and (5 can be x and y. 

We make the observation here, that is discussed in 
more detail in Appendix IB| that the transverse exchange 
interactions can appear in the effective Hamiltonian from 
VCFEs involving only the ground state doublet and first 
excited states because there are nonvanishing , and 
matrix elements between these states. Hence, the ap- 
pearance of these effective transverse exchange interac- 
tions in H^^ is strongly tied to the specific form of the 
Tb2Ti2 07 single ion crystal field wavefunctions. 

Because the bare exchange interaction vanishes if Ii 
and I2 are not nearest neighbors, even if one of the bilin- 
ear interactions in Eq. (|20p is a dipole-dipole interaction. 



only two ions involved in Case C (see Fig. 5(c)), the in- 
teraction falls off as the square of the dipole-dipole inter- 
action: l/jR/j/^l^. 



E. Treatment of the dipole-dipole interactions 

In Section |Vl the semiclassical ground states of the ef- 
fective Hamiltonian arc computed first on a single tetra- 
hedron (Section |VB|), then on a periodic cubic unit cell 
fSection lVPp . In the former case, the bare dipole-dipole 
interaction, i/dd, is truncated beyond nearest neighbors. 
In the latter case, one should not truncate the long-range 
dipole-dipole interaction. This problem has been ap- 
proached in two ways. The first way is to derive the 
effective interactions on a finite but large lattice. Then, 
to obtain the interaction between sublattices a and b on 
a periodic unit cell with sixteen sublattices, the interac- 
tions between sublattice a and all the periodic images of h 
on the large lattice are summed up. A second way to treat 
the long-range dipole is to compute the bare microscopic 
interactions on a periodic cubic unit cell by an Ewald 
summatio n^^i^^ and then to derive the effective Hamilto- 
nian respecting the periodicity. The first approach was 
used in Ref. Here we use the latter. 



F. Summary 

We have now worked out the different types of effective 
interactions that arise in the Hamiltonian i/eff obtained 
from the model of Section FlI Al for Tb2Ti207 to lowest 
order in the virtual crystal field excitations (VCFEs). 
Before describing previously published results obtained 
from the effective Hamiltonian when considering a sin- 
gle (isolated) tetrahedron, we briefly summarize here the 
results of Sections UTT] and IIVI 

The effective Hamiltonian for Tb2Ti207 has been de- 
rived to order ((y)/A) which includes a classical part 
and also interactions coming from VCFEs to lowest order. 
The classical part of the effective Hamiltonian is given by 
VW and is nothing other than the dipolar spin ice model 
(DSIM) with Ising exchange ^classical and dipole-dipole 
couplings that merely differ from those of the microscopic 
model (Section III A[) by a constant factor related to the 
expectation value of the bare angular momentum in the 
crystal field ground states; jZciassicai 

The effective Hamiltonian H^a is expressed in terms of 
spin one-half operators which have different time reversal 
properties (Eq. p4)) ) compared to true angular momen- 
tum operators. A large number of different pseudospin 

(2) 

interactions appear in the quantum term . These 
arc constrained to be invariant under time reversal and 
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to respect lattice symmetries. If we switch off tlie dipole- 
dipole interaction temporarily [V = 0), we find nearest 
neighbor Ising interactions which renormalize those from 
the classical term and, also, transverse terms of the form 
(Tj^fTj^ where a,/3 can be components x or y. Effective 
interactions beyond nearest neighbor are Ising exchange 
interactions between second nearest neighbors and one 
out of the two distinct types of third nearest neighbors 
on the pyrochlore lattice (see Fig. [7]). Finally, three-body 
interactions of the form J^^'^a'a^a^ are also generated. 
When dipole-dipole interactions are restored, {V ^ 0), 
i/off acquires two new types of effective interaction. 

1. New short-range interactions acting between near- 
est neighbors and beyond, decaying as 1/|R|^. 

2. Long range interactions decaying as l/|Rp. 

Both contributions, arising when the dipole-dipole inter- 
actions arc switched on, further renormalize the effec- 
tive nearest neighbor Ising coupling and contribute to 
the transverse effective couplings. 

For later reference we write the nearest neighbor effec- 
tive couplings between ions on sites Ii and I2 as 



ili2 ^1 J2 



(26) 



JZZ -lz ~LZ 

where a and (3 can each equal x and y. 

Now we consider the relative magnitude of some differ- 
ent effective couplings on a lattice which will be of use in 
later sections when we interpret our ground state phase 
diagrams. Fig. [5] shows three different effective Hamilto- 
nian couplings as a function of the bare exchange cou- 
pling jTox when the dipole-dipole coupling T) is fixed at 
the value for TbaTiaOy; V = 0.0315 K. The three cou- 
plings are for nearest neighbor interactions J^^a^a^ and 
jxx~x~x g^j^j ^jjg three-body interaction J^^^tj^o-^o-^, all 
expressed in the local coordinate system. The transverse 
bilinear couplings are averaged over the bonds on a single 
tetrahedron to give an idea of the scale of the interactions 
while the three-body coupling is plotted for bonds con- 
necting sublattices 1, 2 and 3 on a single tetrahedron (see 
Table H]) which is representative of the scale of these in- 
teractions. Looking at these nearest neighbor couplings, 
we see that the Ising interaction changes sign at about 
Jex ~ 0.2 K. For < 0.2 K, the Ising interaction favors 
ice-like order and, for jTox > 0.2 K, it favors all- in/all-out 
ordering (see Fig. 14(a)). The Ising coupling J^" is the 
largest coupling over most of the range of jTox ^ 0.2 K, 
followed by transverse couplings, for instance J^^, and 
then the three-body interaction J^^^. For J'cx ^ 0.2 K, 
the Ising and transverse couplings are of similar magni- 
tude. This is therefore a direct microscopic derivation 
showing the restoration of effective pseudospin isotropy 
that is discussed in Refs. [23 andlM. 



G. Relation with previous results 

In Section IIIII and so far in Section IIVI we have 
presented a derivation of an effective Hamiltonian for 
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FIG. 9: (color online). Plot showing the variation of vari- 
ous couplings in the effective Hamiltonian on a lattice, as a 
function of the bare exchange coupling Jex- The bare Hamil- 
tonian has long-range dipole-dipole interactions. We plot bi- 
linear Ising J"""", transverse J""^ and couplings (Eq. 
between nearest neighbor sublattices 1 and 2 and a three- 
body coupling J"^^ connecting sublattices 1, 2 and 3. The 
difference between the J^^ and J^^ couplings is due to the 
choice of local x and y axes. All couplings refer to pseudospin 
couplings in the local coordinate system. The Ising coupling 
changes sign at about ~ 0.2 K implying a cross-over from 
ice-like order to AIAO order in the absence of other interac- 
tions. For e/ex <; 0.2 K, the Ising and transverse terms are of 
similar magnitude. The three-body coupling is the weakest of 
the three interactions; its variation is shown in the inset. 



Tb2Ti207 to lowest order in the quantum corrections to 
the DSIM - that is, to order {V)/A. If this model is to 
prove useful, it is important to ensure that (V)/ A is not 
so large that higher order terms contribute significantly 
to the low energy physics of Tb2Ti207. To test the asser- 
tion that higher order terms (those of order {{V) /A)" for 
n > 2) are not required, a direct comparison was made 
in Ref. I33! with a microscopic model. This microscopic 
model is the one presented in Section [II Al but with the 
microscopic exchange and dipolar interactions restricted 
to spins on a single tetrahedron and with the crystal field 
spectrum cut off beyond the four lowest energy states. 
More precisely, instead of diagonalizing the full model 
Hcf + V, the single ion crystal field Hamiltonian is diag- 
onalized to obtain states \n) and corresponding energies 
En, whereupon all but the lowest two doublets are ne- 
glected leaving, as the new (truncated, "tr" ) crystal field 
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Hamiltonian. 

^^cf.tr = ^£;o(|l/,)(l/J + |2/J(2,J) 
Ip 

+ {Eo + A)i\3jJ{3jJ + \Aj^){AjJ). (27) 

This model should be a good approximation for sufh- 
cicntly weak interactions given that, when interactions 
are switched on, the excited crystal field levels admix into 
the ground state doublet with the first excited levels hav- 
ing the greatest contribution to the new ground state out 
of all the excited levelsi^^ Including interactions on a sin- 
gle tetrahedron within the basis of four single ion states 
on each site, the model was diagonalizcd exactly for a 
range of bare exchange couplings jTcx- For comparison, 
the effective Hamiltonian was derived on a single tetrahe- 
dron and diagonalizcd computationally to obtain its spec- 
trum. One comparison that has been made from these 
spectra involves looking at the variation in the ground 
state degeneracy as a function of the bare exchange cou- 
pling, jTox, and the crystal field gap as shown in Fig. [31 
As the bare exchange coupling jTcx and the gap A arc 
varied, with V fixed to P = 0.0315 K, there is a phase 
boundary between a region with a singlet ground state 
and a region with a doublet ground state. The boundary 
between these regions is the same for both the effective 
Hamiltonian and the model with a truncated crystal field 
spectrum when the gap A is infinite. The boundaries 
move apart as A decreases. However, the difference be- 
tween the two phase boundaries remains relatively fairly 
small even when the gap is about 18 K, as for Tb2Ti2 07; 
J7gx(A = 18K) for the phase boundaries agree to within 
ten percent. Perhaps most importantly, both calculations 
agree that, as the gap decreases, the region of parameter 
space over which the singlet occurs becomes larger. For 
the parameters estimated for Tb2Ti2 07)^i^ the single 
tetrahedron ground state is a singlet. The ground state 
for the microscopic model, for the Tb2Ti2 07 parameters, 
is mainly a superposition of two-in/two-out states. For 
this reason, this state has been called a quantum spin 

That the singlet ground state of the microscopic model 
shows spin ice-like correlations can be understood from 
the effective Hamiltonian. It is, at first sight, a pecu- 
liar result given the classical DSIM phase diagram de- 
scribed in Section IIIII for which the Tb2Ti207 param- 
eters lie in the AIAO phase. The explanation for this 
behavior lies in the fact that, as the gap A is lowered, 
for fixed Tb2Ti207 bare parameters, the nearest neigh- 
bor Ising exchange coupling J^^ (Eq. ([25]) ) is renormal- 
ized by Ising terms coming from the quantum fluctua- 
tions H^l^ = WnVV (see Section HV]). The variation 
of Ising exchange J^^ due to VCFEs as a function of the 
bare exchange is shown in Fig. llOl when the dipole-dipole 
coupling 2? = 0. The straight line is the part from the 
classical term VVV for which the renormalized (Ising) 
exchange is Jciassicai = JeAJ^) (where (J"") is given in 
Eq. (HI and this formula in derived in Section [111)) . The 



Ising exchange J^^ from the quantum fluctuations is an- 
tiferromagnetic (J^^ > 0) so the sum of the classical 
and quantum Ising couplings is less ferromagnetic (nega- 
tive) than the J^^ without VCFEs. Now consider the ef- 
fect of including nearest neighbor dipole-dipole coupling. 
With the dipole-dipole coupling alone (jTcx = 0), J^^, 
is antiferromagnetic (J^^ > 0), favoring spin ice correla- 
tions (see Section Fill B[) . Indeed, in the Dy2Ti207 and 
Ho2Ti207 spin ice materials, the exchange contribution 
to the Ising J^^ coupling is ferromagnetic but the dipole- 
dipole coupling ensures that the net contribution of the 
bare microscopic couplings to J^^ is antiferromagnetic 
hence frustrating a single tetrahedron and the pyrochlore 
lattice4^ In contrast, we see in the lower part of Fig. [TUl 
that the estimated bare exchange coupling in Tb2Ti207 
of J7ex ~ 0.17K)^ places the classical part of the Ising 
coupling (i.e. the contribution to J^^ from VVV) close 
to zero. But, VCFE corrections lead to an effective an- 
tiferromagnetic coupling J^^ > overall leading to spin 
ice-like correlations for Tb2Ti207. 

In another development, the truncated model with six 
crystal states per site on a tetrahedron was found to ex- 
hibit a magnetization plateau in a [111] field below 100 
mK — further evidence that this model exhibits spin ice- 
like behavior for more antiferromagnetic bare couplings 
than one would expect from the DSIMjS 

In summary, this article is devoted to a derivation of 
a quantum spin-1/2 model with anisotropic interactions 
similar to the models used as starting points in Refs. [3^ 
and H^. From this effective Hamiltonian, we establish 
two things. Firstly, we show that the renormalization of 
the Ising exchange towards antiferromagnetic exchange 
(J^^ > 0) via VCFEs occurs, not only on a single tetra- 
hedron, but also on the full pyrochlore lattice. That is, 
even if, at the classical level (which ignores exited crys- 
tal field levels), Tb2Ti207 could be described by a non- 
frustrated [111] pyrochlore Ising model, virtual crystal 
field excitations can render this system a frustrated spin 
ice one, with additional transverse fluctuations. This is 
the main result of this paper. Secondly, in what follows, 
we study the semiclassical spin correlations that these 
effective interactions produce on a lattice. 



V. SEMICLASSICAL GROUND STATES 

A. Convention and Procedure 

The effective Hamiltonian that was discussed in de- 
tail in the previous sections has a large number of dif- 
ferent effective interactions arising from virtual crystal 
field excitations (VCFEs) from the ground state crys- 
tal field doublet (see Section IIVF[) . As we have seen, 
to lowest (first) order in the perturbation expansion in 
(y)/A, the effective Hamiltonian is an Ising model. The 
lowest order quantum corrections to i/cft include trans- 
verse terms between nearest neighbor spins, three-body 
interactions, and anisotropic interactions extending be- 



16 




...PVP 
■■ PVRVP 
—Total 

0.05 



0.1 



0.15 

(K) 



0.2 



■-■PVP 
■■PVRVP 
— Total 



0.25 




0.25 



0.1 0.15 

Jex (K) 

(b) Isotropic exchange and long-ranged dipole-dipole interaction. 

FIG. 10: (color online). Plot showing how the nearest neigh- 
bor Ising exchange coupling J^^ on a lattice (Eq. (|26|l ') is 
renormalized by the quantum terms of the effective Hamilto- 
nian. The horizontal axis is the bare exchange coupling J^^. 
The top figure is for the Ising coupling when the dipole-dipole 
coupling is set equal to zero and the bottom figure includes 
the dipole-dipole coupling of magnetic ions in Tb2Ti207: 
T) — 0.0315 K. In both figures, the dashed line is the cou- 
pling that appears to lowest order VHV in H^s for a pair of 
neighboring sites. The dot-dash line is the correction obtained 
from the quantum term H^g and the solid line is the sum of 
the two contributions to the Ising exchange. A positive J"^, 
in the absence of other interactions, favors spin ice configura- 
tions on each tetrahedron and a negative sign implies AIAO 
configurations. 



yond nearest neighbors. To gain some understanding of 
the effect of the extra terms generated by VCFEs on the 
physics of this system, it is useful to consider a semi- 
classical spin model with the same interactions as in the 
effective quantum Hamiltonian. 

To obtain the required semiclassical model, we first ob- 
serve that the effective quantum Hamiltonian should be 
written in terms of pseudospins one-half. That is, the ele- 
mentary quantum spins take the form S'" = (1/2)(t" and 
the quantum Hamiltonian couplings derived in the previ- 
ous section are rescaled by a factor of four for the bilinear 
interactions and by a factor of eight for the three-body 
interactions. This model is the most suitable model to 
consider from the point of view of a large spin S expan- 
sion. Once the quantum Hamiltonian is written in terms 
of spins 5", we take these quantum spins into classical 
spins which are vectors of fixed length 5=1/2 with 
components parameterized by spherical polar angles. 

To find the ground states of the semiclassical model, 
we start with a randomly chosen initial spin configura- 
tion on a finite lattice and compute its energy. We then 
make a small random rotation of one of the pseudospins 
and accept this configuration only if the energy of the 
new configuration is lower than that of the old config- 
urations. This procedure is iterated until it converges, 
which happens within 0(10'*) steps. This zero temper- 
ature Monte Carlo may settle into a local rather than a 
global minimum. To alleviate this problem, we repeat 
the process for a number of initial states depending on 
the number of spins treated and look for the minimum 
energy configuration; this also allows us to capture any 
ground state degeneracy. 

Only a small number of independent spins are treated 
in this minimization procedure - we find the ground states 
on a single tetrahedron (four spins) and on a cubic unit 
cell with periodic boundary conditions (sixteen spins). 
With this number of spins, we find that only a small 
number of initial spin configurations O(IO^) is necessary 
in the iteration scheme to find consistency between the 
final energies and to capture discrete degeneracy when 
it arises. However, if there is a continuous degeneracy, 
(as in the XY phase described below), O(IO^) initial spin 
configurations are necessary to confirm its existence. On 
a cubic unit cell, we capture the DSIM ground states 
with ordering vector q = 001 and q = in the limit of 
1/A = 0. When A is finite, VCFEs generate interactions 
beyond the DSIM interactions which may lead to more 
complicated (modulated magnetic moment with incom- 
mensurate wavevector q) ground states might be elimi- 
nated by (an inappropriate choice of) periodic boundary 
conditions. Whether this is indeed the case for the iJoff 
considered below is an open question for future work. 
The key problem we address in computing the ground 
states is to establish whether interactions generated by 
VCFEs beyond nearest neighbor do favor spin ice corre- 
lations over a wider range of jTcx than is observed in the 
absence of such terms. 
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B. Ground states of HcS on a single tetrahedron - 4 
CF states 

The first results that we present are those obtained by 
minimizing the energy on a single tetrahedron to make 
a comparison with the results of exact diagonalization of 
the four crystal field state microscopic quantum model on 
a single tetrahedron described in Section flV Gl The ef- 
fective Hamiltonian required to make the comparison in- 
cludes the nearest neighbor interactions and three-body 
interactions on a single tetrahedron obtained by includ- 
ing only the first excited crystal field doublet in the resol- 
vent (Eq. ([8])) when computing the quantum terms, . 
By truncating all the bare interactions to a single tetrahe- 
dron, the effective Hamiltonian is derived following Sec- 
tion |IV] including all possible ways that the mediating ion 
of Case A in Section HVBI can lie on a single tetrahedron. 
Cases B and C (Sections IIV CI and IIVPP are treated en- 
tirely on the single tetrahedron. At first, we present the 
results when the three-body terms are omitted; we in- 
clude them later. When three-body terms are omitted 
and the derived effective interactions are truncated be- 
yond nearest neighbors, the ground states on a single 
tetrahedron must coincide with those on the full lattice 
under the assumption that the lattice ground states have 
q = ordering wavevector. This is because the inter- 
actions for a q = lattice configuration are the same 
as those on a single tetrahedron except for an extra fac- 
tor of two in the effective couplings on the lattice. This 
factor of two comes from the fact that pairwise effective 
interactions couple a spin on one sublattice, a, to two 
b sublatticcs (a b) whereas, on a tetrahedron, a spin 
with sublattice label a couples to only one b sublattice 
spin. 

The ground states are shown in Fig. 1111 This figure 
shows the energies computed from the effective Hamil- 
tonian with classical spins on a single tetrahedron for 
different (imposed) specific spin configurations and for 
different values of the bare exchange, jTox- The energies 
of the ground states determined by the minimization pro- 
cedure outlined above are shown as well. 

One finds that as the bare exchange coupling jTox (for 
fixed dipole-dipole coupling V = 0.0315 K) becomes 
more antiferromagnetic (i.e. positive and larger), the 
two- in/two-out ground state gives way to an all- in/all- 
out state (Fig. [Ti)) . However, unlike the classical model, 
A = oo, for which only the two-in/two-out and all-in/all- 
out states appear, there is a q = configuration separat- 
ing the two-in/two-out state from the all-in/all-out state 
in which the classical spins lie fully in their local XY 
planes perpendicular to the local [111] directions. 

There is a continuous degeneracy of XY configurations 
such that the vector sum of the spins is zero as one 
should expect for sufficiently strong antiferromagnetic ex- 
change. A single spin configuration among the continu- 
ous set of XY ground states is illustrated in Fig. |14(b)[ 
These ground states belong to the two dimensional irre- 
ducible representation of the point group of the tetrahe- 



dron Oh (see, Refs. l63l andlTOl) which includes the discrete 
ground states of the material Er2Ti207 (see, for exam- 
ple, Ref. [l^ . The onset of the XY phase corresponds 
to a range of values of the bare exchange jTox where the 
effective Ising and transverse couplings (shown in Fig. [9]) 
arc of similar magnitude such that it is energetically fa- 
vorable for the spins to lie in the local XY planes. 

For the classical DSIM with dipole-dipole interactions 
truncated beyond nearest neighbor, there is a phase 
boundary^ between a spin ice state and an all- in/all- 
out phase at jTox = 51? which is roughly 0.158 K for the 
Tb2Ti207 dipolar coupling V = 0.0315 K. Reading from 
the phase diagram in Fig. 1111 the spin ice to XY boundary 
of the nearest neighbor effective Hamiltonian is at about 
0.17 K and the onset of the all-in/all-out (AIAO) phase 
is at about jTox = 0.22 K for the effective Hamiltonian on 
a single tetrahedron. A direct comparison of these num- 
bers with the classical model (DSIM) is possible because 
(i) the phase boundary of the DSIM depends on the ra- 
tio of the effective Ising exchange i/dassicai (Eq- (HH)) to 
the effective dipolc couplings I?ciassicai (Eq- (fTB|) ) and (ii) 
the ratio of effective couplings in the classical term VVV 
is simply the ratio of bare couplings. That the phase 
boundary out of the two-in/two-out Ising configuration 
in the semiclassical effective model appears for more pos- 
itive i/ex than with the classical term {VVV) alone is 
because the effective Ising exchange, J^^, in the quan- 
tum model receives a contribution from H^g that makes 
it more antiferromagnetic (J^^ becomes more negative) 
hence making spin ice Ising configuration energetically 
favorable. We mention also that the effective Ising cou- 
pling, J^^, between Jox = 0.17 K and Jox = 0.20 K is an- 
tiferromagnetic (positive) whereas HcS calculated solely 
to VVV order has ferromagnetic effective Ising exchange 
(favoring all-in/all-out order) over this range. 

When we include the three-body terms on a single 
tetrahedron, we again find three ground state phases with 
the same phase boundaries that we found in the two- 
body case. Whereas the AIAO phase is the same with 
and without three-body terms, the other two two-body 
phases are modified by the introduction of three-body 
couplings. Instead of ground states with spins aligned 
along the local Ising directions, one finds that the spins 
are canted away from the [111] directions while retain- 
ing the spin ice ordering in the Ising components. The 
canting (which is shown in Fig. [T^ for one of the ob- 
served ground states) is such as to preserve the moment 
of the perfectly Ising spin configuration. The variation in 
the canting angle is shown in the inset to Fig. [TT] Only 
in the spin ice regime, J7ex < 0.17 K, are the ground 
states affected by a canting away from the Ising direc- 
tions when three-body interactions are included. The 
three-body terms do not produce a canting away from 
the two-body XY and AIAO ground states. Also, with 
three-body interactions included, the degenerate XY con- 
figurations cease to be the lowest energy states in the in- 
termediate region 0.17 K < J'^x ^ 0.22 K - the continuous 
degeneracy present without three-body terms is broken. 
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The transverse (XY) components of the spins in the spin 
canted state are ordered into six discrete configurations. 
The transverse component configurations, considered on 
their own, are a discrete set of configurations belonging 
to the aforementioned class of XY ground states — they 
are referred to as '02 states in the literaturep^>^'2^ The 
canting angle is zero in the AIAO phase. 

Turning to the ground state energies themselves in the 
presence of three-body interactions, we first note that 
the three-body terms make no contribution to the XY 
and perfectly Ising spin configurations because the three- 
body terms are of the form (t|^ cr"^ (jf^ with a = x,y which 
vanishes in these cases. However, the perfectly Ising-like 
spin ice configurations are not true ground states when 
three-body terms are present for jTox ^ 0.27 K - there is a 
small canting away from the Ising directions. There is a 
difference between the two and three-body ground state 
energies that is small (about 4% at most) arising from 
the relative sizes of the two and three body couplings 
illustrated in Fig. [9l 

Fig. [HI is useful in interpreting the single tetrahedron 
ground state diagrams (with or without three-body inter- 
actions). In both phase diagrams, the sign of the Ising 
coupling determines the correlations of the Ising com- 
ponents of the classical spins (which, as we report in 
Section IV Dl ceases to be true beyond a single tetrahe- 
dron), and the relative magnitude of the Ising coupling 
and other couplings is correlated to the canting of the 
spins away from the Ising directions — the spins being 
furthest from the Ising directions when the transverse 
couplings become comparable to or greater than the Ising 
coupling. 

Now we arc in a position to compare the single tetra- 
hedron results for the quantum four crystal field state 
(ground and first excited doublet) effective Hamiltonian 
with the classical results. The quantum phase dia- 
gram showing the ground state degeneracies is plotted 
in Fig. [3l Focusing on the Tb2Ti207 crystal field gap 
of 1/A = 0.055 K^^, the boundary between the sin- 
glet and doublet states is at about J^-^ = 0.21 K (with 
V = 0.0315 K). The phase boundary for the semiclassi- 
cal ground state derived from the effective Hamiltonian 
(Fig. [TT|) is close to this value, at about jTox = 0.22 K. 
The degeneracies of the semiclassical ground states and 
the quantum ground states on a single tetrahedron agree 
for Jex > 0.22 K whereas for Jox < 0.22 K, the sin- 
glet quantum ground state appears for the same range of 
couplings as both the classical six-fold degenerate canted 
spin ice ground state and the XY phase. 



C. Ground states of HcS on a single tetrahedron - 
13 CF states 

Before leaving the subject of the ground states on a 
tetrahedron, we compute the ground states on a single 
tetrahedron with the full crystal field spectrum included 
in the resolvent operator (rather than considering a trun- 
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FIG. 11: (color online). Semiclassical ground state energy of 
the effective Hamiltonian on a single tetrahedron when quan- 
tum terms are computed including VCFEs only to the first 
excited crystal field doublet (i.e. omitting virtual excitations 
to higher energy crystal field states). In obtaining these re- 
sults, the bare microscopic exchange, -ffcx, and dipolar inter- 
action, -ffdd, were truncated at the nearest neighbor distance. 
The resulting effective couplings generated in were also 
truncated beyond nearest neighbor. The parameters used for 
this calculation are = 0.055 K~\ V = 0.0315 K and 

(J^) = 3.0. Also plotted, are energies of three different im- 
posed spin configurations. The omission of three-body inter- 
actions changes the ground state energy by a few percent de- 
pending on the canting angle produced by these interactions. 
The inset shows the angle from the local z axes of Table |T] 
through which the spins are canted in the ground states of 
the model - for 0.17 K, the canting is due to the three- 

body interactions and the 90 degree canting angles signal the 
onset of the local XY ground states which are ground states 
even without the three-body interactions. 



cation of the spectrum to the ground and first excited 
doublets as we have done in Section FV Bp . In this sub- 
section, the ground states of iJeff we present were com- 
puted as a function of both jTox and A. Here A is an ad- 
justable gap that shifts all the excited crystal field states 
rigidly with respect to the ground state doublet leaving 
the wavefunctions identical to those that one would ob- 
tain by diagonalizing the Tb2Ti207 crystal field Hamil- 
tonian. By artificially varying A in this way, we can tune 
the system from a classical spin ice with nearest neighbor 
dipolar interactions to a model in which VCFEs are sig- 
nificant. The results are shown in Fig. 1131 As one would 
expect based on the limiting case 1/A = of spin ice 
and the results discussed in Section IV Bl for 1/A = 0.055 
K~^, the all-in/all-out ground states and the two-in/two- 
out ground states are separated by a wedge of ordered lo- 
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FIG. 12: The semiclassical LRSIooo ground state of the ef- 
fective Hamiltonian represented on a single tetrahedron. The 
black arrows show the local Ising components of the spins 
and the grey arrows, the (smaller) local XY components. The 
canting does not alter the moment on each tetrahedron rela- 
tive to the moment with uncanted spins. 

cal XY ground states (with continuous degeneracy when 
three-body interactions are omitted). The range of jTox 
over which the wedge extends increases as A decreases. 
For A < 20 K, the AIAO phase is suppressed entirely 
because VCFEs increase transverse couplings relative to 
the Ising couplings (see Fig. (5]). We return to this phase 
diagram in Section IV Dl where we make a comparison of 
Fig. [12] with the ground states on a cubic unit cell with 
periodic boundary conditions. 

D. Ground states with long-range interactions 
included - 13 CF states 

In the foregoing, we have presented the ground states 
for Hcff derived on a single tetrahedron. For the case 
of effective nearest neighbor bilinear spin-spin interac- 
tions, the single tetrahedron ground states are the same 
as the four sublattice (q — 0) ground states on the py- 
rochlore lattice. However, we know that, in the DSIM, 
obtained from the Hcs on a lattice when the crystal field 
gap A is taken to infinity, one of the ground states is a 
sixteen sublattice configuration (with ordering wavevec- 
tor (0, 0,27r/a)) on a conventional cubic unit cell — the 
LRSIooi state shown in Fig. [H The long-ranged nature of 
the dipole-dipole interaction is responsible for the lower 
energy of the LRSIqoi state compared to other ordered 
states that satisfy the local spin ice rules4i This observa- 
tion for the DSIM tells us that we should truncate neither 
the bare dipole interaction to nearest neighbor nor the ef- 
fective interactions and that we should not assume q = 
ordering as was done implicitly in Section IVBI Inspired 
by the case of the DSIM, we investigate the ground states 
on a cubic unit cell with periodic boundary conditions. 
The ground states that we find in this section are for the 
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FIG. 13: (color online). Semiclassical ground states of Hes on 
a periodic cubic unit cell with both bare and effective interac- 
tions truncated beyond nearest neighbor. The phase diagram 
shows the ground states that are obtained over a range of A 
and Jcx- Three body interactions are also neglected. The 
horizontal bar in the top left hand corner represents the un- 
certainty in the bare exchange coupling J'cx in Tb2Ti20 , 
at the value of the crystal field gap A appropriate to this 
material. 




(a) All-in/all-out state. (b) XY configuration. 



FIG. 14: A pair of configurations that are ground states of 
the four sublattice classical analog of the quantum effective 
Hamiltonian without three-body interactions. The all-in/all- 
out state (a) occurs for J^ox > 0.22 K and the XY configura- 
tion (b) for 0.17 < J'cx < 0.22 K. 



pyrochlore lattice, assuming that the magnetic unit cell 
is no bigger than the conventional pyrochlore cubic unit 
cell (with 16 sites). 

For the problem of finding ground states, the effective 
Hamiltonian is derived in the following way, which dif- 
fers from the approach presented above in Sections IVBI 
and IV CI in having to treat the long-ranged dipole-dipole 
interaction. The bare Hamiltonian, which has near- 



20 



LRSI 001 




^■^0 0.05 0.1 0.15 0.2 0.25 

Jex (K) 

FIG. 15: (color online). Ground state energy for an effective 
Hamiltonian derived from a model with isotropic exchange 
and long ranged dipole-dipole interactions on a cubic unit cell 
with periodic boundary conditions treated by an Ewald sum- 
mation with = 0.055 K"^ and V = 0.0315 K. The ground 
states with (circles) and without (squares) including three- 
body interactions are shown as well as the energies of different 
(imposed) spin configurations. The two-body ground states 
are the LRSIooi state for ^ 0.06 K which is a two-in/two- 
out state with ordering wavevector (0, 0, 27r/a) and an or- 
dered two-in/two-out state with ordering wavevector (0,0,0) 
for -SJ 0.06 K. When three-body terms are included the 
spins cant out of the Ising directions as indicated in the inset. 



est neighbor isotropic exchange and long-ranged dipole- 
dipole interactions, is computed on a sixteen site cubic 
unit cell by summing the dipole-dipole interaction over 
all periodic images by an Ewald summation.— The ef- 
fective Hamiltonian is then computed numerically for this 
periodic model on a cubic unit cell (using identities at the 
end of Appendix |B] and summing over the full 13 state 
crystal field spectrum in the perturbation theory). This 
procedure preserves the periodicity of the Hamiltonian. 
One could have instead derived the effective Hamiltonian 
on a lattice and then sum the interactions over a large 
but finite lattice assuming periodicity in the classical spin 
configurations on a cubic unit cell. We present results, in 
this section, for the former case, but the latter approach 
gives results that are quantitatively very similar.— 

The semiclassical ground states energies of the result- 
ing model are computed by replacing the pseudospin op- 
erators (1/2)0-" with classical spin components S*". They 
are given in Fig. [T3] as the bare exchange coupling is 
varied with 1/A = 0.055 K"! and V = 0.0315 K; the 
values appropriate to Tb2Ti207. In the same figure, we 
also plot, for comparison, the energies of various imposed 
(fixed) spin configurations. As in the single tetrahedron 



case, it is useful to distinguish the ground states ob- 
tained when three-body spin interactions are removed 
and the ground states for the model with all interactions 
included. 

Without three-body spin interactions (open squares in 
the main panel of Fig. [T5|) we find that for weakly an- 
tiferromagnetic J'^^ < 0.06 K, the ground state is the 
LRSIooi state. But, for more antiferromagnetic jZoxj in- 
stead of the all-in/all-out state found for the DSIM (see 
inset of Fig. [2]), the ground state, at least up to Jox = 0.4 
K, is a state with identically ordered tetrahedra (order- 
ing wavevector q = 0) each obeying the ice rules with 
spins in the local Ising directions — we refer to this as 
the LRSIqoo phase. 

If the dipole-dipole interaction is cut off at nearest 
neighbor in the microscopic bare model of Section III Al 
before computing the effective Hamiltonian and if effec- 
tive couplings beyond nearest neighbors are removed, the 
LRSIooi Ising state has the same energy as the LRSIooo 
Ising phase for jTex ^ 0.06 K. The XY phase that we 
found on a single tetrahedron does not appear in the 
conventional cubic unit cell model unless all effective in- 
teractions are cut off beyond nearest neighbors. 

When three-body interactions are restored (filled cir- 
cles in Fig. I15p , the spins cant away from the Ising direc- 
tions and the energies are lowered relative to the ground 
states with only two-body interactions considered, sim- 
ilarly to what was found in Section IV Bl on a tetrahe- 
dron. The ordering of the Ising components of the spins 
is the same regardless of whether three-body interactions 
are present or not. The canting angle of the spins away 
from the Ising directions is shown as the inset in Fig. [151 
There is a maximum in the angle at about J^ex = 0.11 
K and two minima at about 0.02 K and 0.21 K over 
the explored range of jTcx- The greatest angle is about 
8 degrees, compared to 90 degrees in the case of a sin- 
gle tetrahedron. Broadly, the variation in the canting 
angle follows the magnitude of the three-body coupling 
(shown for a choice of three sublattices in Fig. [5]). The 
minimum in the canting angle at about J^^x ^ 0.02 K 
coincides with a minimum in the mean squared three- 
body coupling over all such couplings on the cubic unit 
cell at this value of the bare exchange. The minimum 
implies that there is little energy gain to a canting of 
the spins. The non-monotonic change in the three-body 
couplings is allowed because the coupling has contribu- 
tions both quadratic and linear in jTox coinciding with the 
exchange-exchange and exchange-dipole contributions to 
HoS of Eq. (|20p . The second minimum in the canting 
angle, jTox ~ 0.21 K, coincides roughly with a vanishing 
in the three-body coupling (shown in the inset to Fig. [9]) 
and with a change in the sign of the nearest neighbor 
Ising coupling. We expect therefore, two effects to be 
at work - a weakening of the three-body coupling and 
the same effect that suppressed the canting angle in the 
transition from the XY phase to the AIAO phase on a 
single tetrahedron (Section IV Bp . The difference in this 
case is that it is a weak effect compared to that of the 
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effective further neighbor interactions which control the 
Ising ordering in this range of coupHngs. The ordering of 
the local XY components of the spins is identical to that 
described in Section IVBI - one such type of XY ordering 
is shown in Fig. [12] in which smaller arrows indicate the 
canting direction away from the Ising directions.— 

We note in passing that the material Tb2Sn207 which 
is, microscopically, very similar to Tb2Ti207>^ under- 
goes a phase transition at about 0.87 K to a magnet- 
ically long-range ordered phase with ordering wavevcc- 
tor q = and two- in/two-out spin ice configurations on 
each tetrahedron 1^ The spins in this ordered phase are 
canted away from the local Ising directions. In these 
respects, the spin canted LRSIqoo ground state of the 
effective Hamiltonian discussed above is similar to the 
magnetic order in Tb2Sn2 07. But the nature of the spin 
canting differs between the model and the material. The 
spin canting in Tb2Sn2 07 is such as to reduce the mo- 
ment, on each tetrahedron, compared to the moment if 
the spins were not canted^^ whereas the canting of the ef- 
fective spin 1/2 in the LRSIqoo state (indicated in Fig.fT^ 
produced by the three-body terms gives a moment, on 
each tetrahedron, that is the same as the moment of the 
LRSIooo configuration without spin canting. 

The classical DSIM (which is recovered for an infinite 
ground to first excited crystal field gap, A) with long 
range dipole-dipole interactions has a phase boundary be- 
tween the LRSIooi spin ice configurations and all-in/all- 
out states at about J^ex = 0.14 K (see inset of Fig. [2]). We 
have seen that, in the effective model of Tb2Ti207, the 
semiclassical ground states, at least for J'ex < 0.25 K, 
are spin ice configurations although the bare exchange 
coupling J7ex is antiferromagnetic so we see that spin 
ice correlations are favored by VCFEs. However, Fig. [9] 
shows that the average nearest neighbor Ising exchange 
J^^ swaps sign at about 0.2 K so the persistence of ice- 
like correlations in Tb2Ti207, in the form of the LRSIooo 
state, up to, at least, J7cx = 0.4 K (see Fig. [2]) is not due 
to the renormalization of the Ising exchange described 
in Section IIV Gl but is induced by further neighbor cou- 
plings. The effective Hamiltonian to order A(j7cx/A)^ 
is therefore a novel two-in/two-out model that does not 
rely on nearest neighbor interactions to produce ice-like 
correlations. 

To shed some light on the fact that the all- in/ all- out 
state, (Fig. 14(a) ), observed on a single tetrahedron and 
in the DSIMj^i^ii^ is not seen in the sixteen sublat- 
tice case, (for the value 1/A w 0.055 K"-'^ as shown in 



Fig. [2] and Fig. [T5)) . we have computed the semiclassical 
ground states for a range of A and bare exchange cou- 
plings for the sixteen sublattice effective Hamiltonian on 
a cubic unit cell with periodic boundary conditions. We 
have omitted the three spin interactions which are not 
responsible for the presence of the LRSIqoo state. The 
result is shown in Fig. [21 This is the main result of our 
paper. 

For A > 340 K (1/A < 0.003), the phases arc those 
of the DSIM with a phase boundary at about jTox — 



0.14 K when 1/A = 0. For comparison, we include an 
inset showing the classical DSIM phase diagram for V = 
0.0315 K4iii^ As the gap A is lowered from infinity, a 
q = LRSIooo phase — appears at about A ~ 340 Kj^ 
As A is lowered further, the range of jTox over which 
this LRSIqoo phase is observed increases — the spin-spin 
interactions arising from VCFEs stabilizing the LRSIqoo 
state. Indeed, for A < 29 K (1/A - 0.035 K^^ ), there 
is no all- in/all-out phase at least for any jTox < 0.4 K. 

Over the range of jTox explored here, the LRSIqqq phase 
boundary has a dip with a minimum at about jTox = 0.14 
K. The LRSIqoo is not observed in the DSIM so the quan- 

(2) 

turn terms of H^^f is responsible for its existence. There- 
fore J'cx = 0.14 K is the exchange coupling at which 
the effect of the classical term is minimized because the 
isotropic exchange and the dipole-dipole contributions 
to the Ising exchange almost cancel each other. This 
accounts for the "tail" in Fig [2| where the LRSIqqi to 
LRSIqoo phase boundary extends to 1/A ^ 0.005 — 
the quantum terms are dominant at about jTox = 0.14 K. 
The phase diagram in Fig. [2| shows that for jTcx 0.25 
K, the LRSIqoo spin ice appears over a larger range of 
1/A. This is because for larger values of the exchange, 
the quantum terms are larger for a given A and also 
because the quantum terms vary as they eventually 
dominate over the classical terms. 

These observations lead us to two comments. Firstly, 
the shape of the LRSIqqq phase boundary in Fig. [51 is 
similar to the shape of the phase boundary for nearest 
neighbor bare and effective interactions shown in Fig. 1131 
The explanation we have given earlier in this Section FVDI 
(for the case with long-range dipoles on a cubic unit cell) 
for the shape of this boundary is equally applicable to 
the case with nearest neighbor bare and effective inter- 
actions discussed in Section IV CI A comparison of these 
two figures reveals that, whereas quantum terms strongly 
influence the nearest neighbor phase diagram enough to 
produce an XY phase, the replacement of this phase by 
LRSIqoo requires interactions beyond nearest neighbor 
which, therefore, should not be neglected. 

Secondly, as J7cx/A increases, eventually higher order 
terms in powers of J7cx/A must become important and 
our effective model will break down. It is possible that 
the inclusion of higher order terms would lead to the all- 
in/all-out phase persisting to larger values of 1/A than 
we find considering only the lowest order quantum cor- 
rections -ff^^' to the DSIM. Fig. [3| is a comparison of 
the exact four state model of Eq. (|771) with the effective 
Hamiltonian on a single tetrahedron. The singlet-doublet 
phase boundary indicates that the AIAO phase region 
should occupy a larger range of 1/A as jTox increases 
than is borne out by the semiclassical ground state cal- 
culation on a single tetrahedron (Fig. [TT|) . On the basis 
of this comparison alone, however, one cannot draw any 
conclusions about the effect of higher order corrections 
on the phase diagram of the effective Hamiltonian on a 
lattice. 
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VI. SUMMARY AND DISCUSSION 

In this paper, wc have introduced a low energy effec- 
tive Hamiltonian for Tb2Ti207 formaUy derived from a 
minimal microscopic bare Hamiltonian. The bare model 
consists of a crystal field Hamiltonian for each magnetic 
ion and isotropic exchange and dipole-dipole interactions 
between the J = 6 angular momenta (Eq. ([1]) and discus- 
sion m Section HT^ . The low energy model is expressed 
in terms of effective spin one- half operators which operate 
on states in the two dimensional Hilbert space spanned 
by the single ion ground state crystal field doublet on 
each magnetic site. The effective theory is obtained as a 
perturbation expansion in {V) / IS. where {V) is the char- 
acteristic energy scale of the spin-spin interactions which 
incorporate exchange and dipole-dipole coupling;^ and 
A is the energy gap between the ground and first excited 
levels in the crystal field spectrum. In Section[TTl we gave 
a detailed discussion of the terms that arise in the effec- 
tive model to lowest order in the quantum corrections. 
To first order in the effective Hamiltonian in powers of 
(F)/A is the DSIMiiii^i^ which has only interactions 
between the Ising components of the pseudospins. This 
model on its own has an antifcrromagnetic all- in/all-out 
(AIAO) ground state for the estimated bare exchange, 
dipolar and crystal field parameters for Tb2Ti207 (see 
vertical dashed line in the inset to Fig. [2|). The next 
(second) order in {V) / IS. includes the lowest order quan- 
tum fluctuations involving virtual transitions into excited 
crystal field levels. We found that the introduction of 
these virtual fluctuations leads to a renormalization of 
the effective Ising exchange coupling in i?cff of the lowest 
order (spin ice) model in such a way that two-in/two- 
out Ising configurations are favored on the single tetra- 
hedron over a wider range of jTox than one would find 
from the lowest order (DSIM) model. Also, to second 
order in (y)/A, various anisotropic transverse effective 
exchange couplings appear (in addition to corrections to 
the effective Ising exchange) as well as some three-body 
interactions. Broadly speaking, the interactions between 
the effective spins become less Ising-like in the presence 
of virtual crystal field excitations (VCFEs). This be- 
havior is also borne out by comparisons of the diffuse 
neutron scattering pattern for Tb2Ti20 , ^^iS?, 28^53 
both classical mean field theory with classical Heisenberg 
spins and finite Ising-like anisotropy^ and by RPA cal- 
culations starting from the bare microscopic model pre- 
sented in Secti on III A\r ^ In other words, the conclusion 
reached in Ref. l33l that, on the basis of exact diagonal- 
ization calculations and perturbation theory calculations 
on a single tetrahedron, Tb2Ti2 07 may perhaps be de- 
scribed by a soft (quantum) spin ice system is upheld by 
the work presented in the present paper. 

We studied the properties of the low energy (effec- 
tive) Hamiltonian ifeff by finding the ground states, as a 
function of bare isotropic exchange couplings J^ex (from 
the model in Eq. Q), and for the crystal field spec- 
trum of Tb2Ti207, under the assumption that the pseu- 



dospins are classical (i.e. the pseudospins are vectors 
of fixed length S = 1/2). Truncating the bare Hamilto- 
nian and then the effective Hamiltonian to nearest neigh- 
bor interactions and assuming ground states with order- 
ing wavevector of q = (identical spin configurations 
on elementary tetrahedra on the pyrochlore lattice) and 
omitting three spin interactions, we found three different 
scmiclassical ground states depending on the ratio J^o^/'D 
(see Eq. (HSl)). Specifically, we found (i) a two- in/two-out 
state and (ii) an all- in/all-out state. In addition to these 
two states is one with spins lying in the local XY planes 
perpendicular to the [111] directions (see, for example, 
Fig. |14(b)P with a continuous degeneracy. 

If, instead, the original model has long-ranged dipolar 
interactions treated by an Ewald summation on a single 
cubic unit cell, then the effective Hamiltonian has inter- 
actions with the periodicity of a cubic unit cell. For such 
a model, again without three spin interactions, (and as- 
suming that the magnetic unit cell is not larger than a 
single cubic unit cell), wc find (for 1/A ~ 0.055 K"-'^ rel- 
evant to Tb2Ti207) two scmiclassical ground states. For 
weakly antifcrromagnetic bare exchange, j/ox , the ground 
state is the LRSIqoi phase (see Fig.[T]) — a ground state of 
the dipolar ice model — and, for more antifcrromagnetic 
J7cx, the ground state is a two-in/two-out state with prop- 
agation (ordering) wavevector q = 0. The latter result 
— the persistence of spin ice correlations with antifcrro- 
magnetic bare coupling — is partly a consequence of the 
renormalization of the effective Ising exchange coupling 
which includes contributions from the bare dipole cou- 
pling T) and the bare isotropic exchange coupling jTox- 
It is also partly due to the presence of further neighbor 
interactions not present in the microscopic model. Be- 
cause spin ice-like correlations appear over a wider range 
of couplings than one would find in the classical model, 
the VCFEs are responsible for frustrating the interac- 
tions in our simplified model (see Eqs. and ([5])) of 
Tb2Ti207. 

When the three-body interactions are incorporated, 
the ordering of the Ising components of the spins is not 
changed from the results without three-body terms, but 
the effective spins then become canted out of the local 
Ising directions and the local XY components are ordered 
into the so-called -02 states (see Fig. This XY 

ordering is observed in the easy plane antifcrromagnetic 
Er2Ti207. However, we note that the effective Hamilto- 
nian for Er2Ti2 07 has no three-body interactions (a con- 
sequence of time reversal within a Kramers doublet) so 
the effective Hamiltonian for Er2Ti2 07 cannot account 
for the observed ordered state by means of three-body 
interactions. 

Returning to Tb2Ti207, in the present work we have 
established that VCFEs can be included as a significant 
perturbation to the DSIM and that the interactions in- 
duced by VCFEs have an important effect on the physics 
of this material. By far the most important problems 
now remaining are to establish the ground state and low 
energy excitations of the fully quantum effective Hamil- 
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tonian derived in this paper beyond the single tetrahe- 
dron approximation (see Section HVGp and to assess the 
importance of higher order terms in the perturbation ex- 
pansion. This might be accomphshed by pursuing exact 
diagonahzation or series expansionsi^ 

Further unresolved problems are to account for the 
long-range order in Tb2Ti207 that is induced by a [110] 
magnetic ficld^ and by applying pressure There is also 
evidence to suggest that there are dynamical lattice dis- 
tortions away from a pyrochlore structure in zero mag- 
netic field^i^ However, the extent to which these af- 
fect or are affected by the magnetism in the material is 
not known. With the availability of an effective Hamilto- 
nian that considers the effect of excited crystal field levels 
in Tb2Ti207, one can perhaps hope to supplement the 
model to explore the role of the lattice on the magnetism 
of Tb2Ti207 and Tb2Sn207. 

One could extend the work in this article by including 
interactions in the bare Hamiltonian besides isotropic ex- 
change and dipole-dipole interactions. For example, one 
could explore the effect of generalized anisotropic nearest 
neighbor exchange interactions as was done at the mean 
field level in Ref. M (for Yb2Ti207) and Ref. M (for 
Er2Ti2 07). In addition, one could include further neigh- 
bor interactions in the bare Hamiltonian. It is already 
known that further neighbor interactions are present in 
the related (spin ice) material Dy2Ti2074i If further 
neighbor interactions were shown to be significant in 
Tb2Ti2 07, they could be incorporated following the ap- 
proach in this article. 

Looking beyond the question of the ground state of 
Tb2Ti2 07, we point out that an effective Hamiltonian 
of the type described in this article might find a use in 
other problems on magnetic systems. For example, this 
approach might find some use in studying the material 
Pr2Sn2 07 -SS. which has been referred to as "dynamic 
spin ice" with an ill-understood fast dynamics compared 
to Ho2Ti207. Two other pyrochlore magnets with Ising- 
like crystal field ground states are the metallic spin ice 
Pr2lr207 — and the material Pr2Zr207— both of which 
exhibit no long range magnetic order at low tempera- 
ture. Finally, we mention another material for which the 
effective Hamiltonian formalism might be useful — the 
langasites Nd3Ga5SiOi4^'2ii92^ ^^^^ PrgGagSiOM--^'^^ 
These materials show no sign of order at least down to 35 
mK although the scale of the interactions in both com- 
pounds is much larger, as read off from the Curie- Weiss 
temperatures (-52 K and -2.3 K for Nd3Ga5SiOi4^i22 
and PraGasSiOi^S^ respectively). 

We hope that the present work stimulates further the- 
oretical investigation into the exotic and interesting be- 
havior displayed by these materials. 
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APPENDIX A: EFFECTIVE HAMILTONIAN 

In order to keep this paper self-contained, and in the 
hope that the approach we have followed here will be of 
use to others, we sketch out the main ideas behind the 
derivation of the effective Hamiltonian. The discussion, 
which we keep fairly general, roughly follows Ref. [6^ to 
which we refer for a broader context. 

We consider a quantum mechanical system described 
by Hamiltonian H which can be split into Hq plus a small 
perturbation V. We label the exact eigenstates of H by 
l^'n) which corresponds to eigenvalue £"„ for n from 1 to 
the dimension of the Hilbert space Af. The eigenstates of 
the Hamiltonian Hq are denoted jrio) (where the integer 
n labels different eigenstates) and satisfy 

HoIuq) = So,„|no). 

In the following, we imagine that the ground state of 
Hq is p-fold degenerate and that eigenstates are labeled 
|lo) to Ipo) and have energy Eq. When we introduce the 
perturbation V, to zeroth order in ordinary degenerate 
perturbation theory, the ground state wavefunctions are 
some particular admixtures of these degenerate states — 
in this sense they are strongly coupled by the perturba- 
tion. 

We wish to set up a Hamiltonian that "lives" in the 
subspace spanned by the ground state levels of Hq and 
which includes the effect of V on these levels. There- 
fore, we introduce a projector V that projects onto this 
subspace. Given an exact state 

T'l^'™) = |*0,n) 

where I'I'cti) is a linear combination of |no) for n = 
1, . . . ,p. We refer to this subspace as the model space 
Because the perturbation is assumed to be weak, the 
exact eigenstates |\E'n), for n from 1 to p lie mainly within 
DJl. We also introduce an operator fl that "undoes" the 
effect of the projector V, 

It follows that 1^*0, n) = T'f^l^o.Ti) and, because this equa- 
tion is satisfied by any linear combination of the exact 
states |*o,«), we find that V^V = V. 
The following intermediate result holds: 

[VL, Ho]v = vnv - nvvfiv. (Ai) 

To see this, begin with the Schrodingcr equation in the 
form {En - Ho)\^n) = and multiply (from the 
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left) by V to get 

because the projector commutes with the Hamihonian 
The effective Hamiltonian, i?off, is defined to be 

Hcff = VHnV (A2) 
which has the property 

-ffcff|*0,n> = En\^o,n)- 

The effective Hamihonian has eigenstates hving in the 
model space 9Jl and has as eigenvalues the exact eigen- 
values. The projector on the right-hand-side is there to 
ensure that the remaining operators VHQ. operate on 
the model space, 9JT. Operationally, in Eq. (|A2[) . the 
Q. operator rotates the model space state into an exact 
eigenstate. H produces the exact eigenvalue and then the 
exact eigenstate is projected back into the model space. 

We compute iJcff in perturbation theory by expanding 
51 implicitly in powers of V 

O = f + +^(2) + , ^ (A3) 

It is then possible to eliminate O^*"') by introducing the 
so-called resolvent operator 

where Q = 1 — V . The resolvent has the spectral repre- 
sentation 

Note that 7^ = 7^Q. 

To eliminate introduce the series (jA3[) into iden- 
tity (|Aip and use the fact that V projects onto states 
with the same Hq eigenvalue Eq to obtain 

= -JIVV 

and so on. These recursion relations can be solved to get 
Q.^^'i in terms of 7^, F and "P. 

The effective Hamiltonian ()A2|) then takes the form 

iJcff = VHV + VHUHV + ... 
which, in turn, is 

H^ff = vHoV + vvv + vvnvv + . . . 

because, in term VHTZHV, the unperturbed Hamilto- 
nian is eliminated because it does not contain any terms 
that connect the model space with the space orthogonal 
to it — that is, terms like VRqIZHV vanish. 



APPENDIX B: CALCULATIONS FOR CASE A 

In this appendix, we give more details of the calcula- 
tion leading to the effective pseudospin interactions from 
i/^ff which is that part of the effective Hamiltonian that 
includes VCFEs to lowest order in {V) /A. We begin with 
Eq. pTj) which we reproduce below 

E EE^("^i'"^2,m3)^f,,Jr,j£ 

Q,/3,p,(T nip W 

^ \mi^ ,4, Wi^ , TOj3,3) ("^Jl ,4, Wi^ , TO/3,3 1 
Eq — Ew 
X iCPl.JlJlVimi, m,,me). (BI) 

In this formula, the lattice sites /i, I2 and It, have been 
fixed. We observe that the matrix elements for the angu- 
lar momenta on sites Ii and /3 arc taken between states 
within the ground state crystal field doublet. The nonva- 
nishing matrix elements within this doublet are given in 
Eq. (j4|). From this equation, we see that a and a must 
equal z. We consider the operators in Eq. (jBip that act 
on site Ii 

E |toi)(toi| JfJm4)(TO4| 

mi ,m4 

— >|1)(1|(1|J;J1) + |2)(2|(2|J|J2) ^ 

= -(.7fJ(|l)(l|-|2)(2|)^-(J/Ja,V 

A similar calculation gives —{J^)a^ on site 13. Eq. (jBip 
becomes 

( ^ M ,("^2|jfjVl^)(l^|JfJm5)\ 

X E E • 

(B2) 

The sum over W runs over all single ion crystal field 
excited states. We will consider only the sum over the 
lowest excited crystal field doublet states: |3) and |4). 
The denominator E\y — Eq equals A. The relevant matrix 
elements are, from exact diagonalization of the crystal 
field Hamiltonian, 



(1|J-|3) 


= A 


(1|>|3) 


= ~iA 


(1|J^|4) 


EE B 


(2|J13) 


= -B 


(2|jn4) 


= -A 


{2\jy\A) 


EE -iA. 



All other matrix elements vanish. We shall not make any 
assumptions about the form of K."^ except for symme- 
try under swapping both pairs of indices. The sums in 
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brackets in Eq. (IB2p give the operators 

ii)(2| (-^ir/,^i,% +^^?f,,^&3 -^fr/,^?^,3 
+i2)(i| - - ^r^iMi. 

((^ir/.^?.\ +^?r//&3 +*^fr/.^f:/3 
+i2)(2i ((^ir,,^f,\ +^jf,/?:,3 -*^ir/.^/:/3 

+z^^f,/f,\,) . (B3) 

Of these four operators, the top two involve virtual exci- 
tations on ion I2 that do not return the ion to its orig- 
inal state but instead take it into the other crystal field 
ground state on ion I2 — an overall Ising spin flip. We 
shall see that, as we should expect, these spin flip oper- 
ations correspond to or effective operators. This 
leads us to an important point - in order for cr^ or d'^ ef- 
fective operators to be significant in the effective Hamil- 
tonian for Tb2Ti207, there must be nonvanishing J^, 
and J^' matrix elements between the ground state dou- 
blet and first excited doublet. In order for this to be the 
case, the ground state and first excited wavefunctions, 
which have the form 

\n)= J2 l^'^-^> 

M=-J 

cannot have only the predominant | J, ±4) (in the ground 
doublet) and |J, ±5) (in the first excited doublet) coef- 
ficients for then the matrix elements would vanish. 
Hence the conclusions of the paper are unlikely to carry 
over to other materials, for example, to the spin ices. 

Referring to Eq. (|f ip . we find that the above opera- 
tors in Eq. (|B3[) can be re-expressed in terms of Pauli 
matrices. So the result of the sum of excited crystal field 
states in Eq. (|B2p is 

-{'<^ZAs+!^KiMis)^Bai- 

(B4) 

After substituting the couplings K. we find that the tr^ 
terms vanish. The resulting expression is time-reversal 
invariant. Incorporating Eq. (jB4p into Eq. (jB2p . we 
find that the overall interactions are, as we stated in the 
main text, Ising interactions cr|^ af_^ (arising from the unit 
operator in Eq. jEH) and three-body interactions of 



the form aj_^af^af^ with a = x,y. Having determined 
the general form of the interactions and their couplings, 
we carry out a sum over all lattice sites Ji, I2 and 73. 
The calculations for Cases B and C in the main text are 
carried out in a similar manner. 

In order to organize the calculation of the terms in 
the effective Hamiltonian, all sums over virtual excited 
states and lattice sites are carried out numerically and 
the calculations described in this appendix are performed 
by exploiting the orthogonality of the Pauli matrices. As 
an example, suppose that the operator coefficients in Eq. 
(|Bip have been evaluated in the |I), |2) basis. We call 
this operator O. We want to decompose this operator 
into a sum of the form 

a,b,c 

where the sum runs over the Pauli operators ct^, ct^, 
(T^,and the unit operator. Coefficients Aabc are deter- 
mined from 

o 

This formula is sufficient for Case A (Section IIV Bp with 
operators on three pyrochlore sites Ji, I2 and 1^. For 
Cases B and C (Sections IIV CI and IIV Pp . we decompose 
into a sum 

a,b 

using 

Bab = \Tr[da''a\ 

APPENDIX C: CRYSTAL FIELD PARAMETERS 
FOR Tb2Ti207 

The crystal field parameters for Tb2Ti207 are obtained 
from those found for Ho2Ti207 in Ref. \^ from the for- 
mula [3] Ref. [43 uses the convention 

' / /I \ 1/2 

for the crystal field parameters. One can convert the set 
of to the using the parameters given in Ref. [o^ 
and the matrix elements of Table [III 

The Ho2Ti2 07 crystal field parameters are: 

7^22- = 791K 7^^ = 3189K 7^^ = I007K 

(i>2)Ho (04)Ho (06)Ho 

7^^ = 739K ^ -725K 7^^^ ^ 1179K. 

i>-'4)Ho (OejHo (OejHo 

(CI) 

The radial expectation values (r™) are given in Table ITTl 
— and the Stevens factors for Tb2Ti207 are given in 
Table HIP 
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TABLE II: Table of radial expectation values, (r™), for 
TbaTiaOy and HoaTizOr.— 

(r^) (r^) 

Ho 0.7446 1.3790 5.3790 

Tb 0.8220 1.6510 6.8520 



TABLE III: 


Table of 


Stevens factors for 


TbaTiaOr and 


HoaTiaOr.-^ 








R3+ 


S2(xl02) 


S4(X10^) 


S'6(X10«) 


Ho 


-0.2222 


-0.3330 


-1.2937 


Tb 


-1.0101 


1.2244 


-1.1212 
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